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1.  Introduction 


During  the  grant  period,  the  research  focus  was  to  use  existing  computational 
tools  and  explore  atomic-scale  dynamics,  defect  reactions,  and  diffusion  in  solids. 
These  tools  include  state-of-the-art  first-principles  calculations  (density  functional 
theory  in  the  local  density  approximation  for  exchange-correlation  with 
pseudopotentials  and  supercells)  as  well  as  classical  molecular  dynamics. 

in  particular,  we  addressed:  a)  the  dynamics  of  dopant  impurities  in  Si  in  the 
heavily-doped  limit,  including  dopant  deactivation  through  clustering,  thermal 
reactivation,  and  enhanced  diffusion;  b)  the  segregation  of  dopants  in  Si  grain 
boundaries  and  dislocation,  where  we  predicted  selective  segregation  in  particular 
columns  in  the  form  of  pairs  or  ordered  chains  ~  the  prediction  was  subsequently 
verified  by  atomic-resolution  transmission  electron  microscopy:  c)  the  mechanism  of 
oxygen  migration  in  crystalline  Si  where  we  found  enormous  complexity  and 
coupled  barriers  in  a  seemingly  very  simple  atomic  jump;  d)  the  flow  of  heat  across 
a  grain  boundary  using  molecular  dynamics  which  led  to  a  direct  calculation  of  the 
Kapitza  resistance  of  a  grain  boundary;  and  e)  mesoscopic  calculations  of  void 
grov^rth  under  electromigration  conditions. 

Brief  descriptions  of  the  main  results  are  given  below.  More  detailed  discussion 
is  given  in  the  attached  Appendices. 

2.  Dynamical  phenomena  in  heavily  As  doped  Si 

In  the  dilute  limit.  As  impurities  occupy  substitutional  sites  and  act  as  donors.  It 
has  been  possible  by  non-equilibrium  techniques  to  incorporate  As  impurities  as 
substitutional  donors  at  very  high  concentrations,  as  high  as  10^^  cm‘^.  However, 
moderate  heating  (~300-400°C)  leads  to  rapid  deactivation.  Partial  reactivation  is 
possible  by  heating  to  higher  temperatures  (>800°C).  Early  work  attributed  the 
deactivation  to  the  formation  of  ASnVm  clusters  (where  V  stands  for  a  vacancy),  but 
there  has  been  no  consensus  about  the  values  of  m  and  n  for  the  dominant  defect. 
Another  major  puzzle  has  been  the  fact  that,  in  the  dilute  limit,  there  is  no 
observable  diffusion  in  the  300-400°C  temperature  range.  Thus,  some  other 
diffusion  mechanism  must  be  operating  in  the  heavily-doped  limit  to  enable  the  As 
atoms  to  migrate  and  form  clusters.  Early  diffusion  measurements  in  heavily-doped 
samples,  however,  found  that  the  diffusion  constant  increases  gradually  with 
increasing  As  concentration,  but  then  begins  to  drop  for  As  concentrations  above 
~3x10^°  cm'^.  More  recent  experiments  using  rapid  thermal  annealing  techniques. 


found  a  sharp  increase  in  the  value  of  the  diffusion  constant  beginning  at  the  same 
concentration.  References  to  all  the  relevant  experimental  work  can  be  found  in 
Appendix  A.  Though  explanations  for  some  of  the  observations  have  been 
proposed  by  earlier  workers,  no  comprehensive  account  of  all  the  data  has  been 
possible. 

A  set  of  comprehensive  calculations  carried  out  by  Ramamoorthy  and 
Pantelides  provided  a  systematic  account  of  all  the  perplexing  observations  and 
made  additional  predictions.  The  detailed  results  are  discussed  in  Appendix  A.  In 
particular,  it  was  found  that: 

i)  Complexes  of  the  form  AsnVm  have  a  negative  formation  energy  for  n>2.  For  all 
such  complexes,  the  binding  energy  of  each  As  atom  (relative  to  an  isolated 
substitutional  site  in  the  bulk  crystal)  is  -1.5  eV.  These  results  mean  that  As- 
doped  Si  is  a  metastable  system  at  all  As  concentrations.  If  it  were  not  for 
kinetic  constraints,  all  As  atoms  would  decorate  vacancies  and  multivacancies. 
The  result  is  a  consequence  of  the  chemical  valency  of  As  which  makes 
threefold-coordinated  As  energetically  favored. 

ii)  It  has  long  been  known  that,  in  the  dilute  limit,  As  diffusion  is  mediated  by  AsV 
pairs.  The  pair  migrates  as  follows:  The  vacancy  loops  around  a  six-member 
ring  and  positions  itself  “in  front”  of  the  As;  the  As  atom  jumps  into  the  vacant 
site;  the  vacancy  again  loops  around  another  six-member  ring  and  again 
positions  itself  “in  front”  of  the  As  atom.  The  activation  energy  for  diffusion  by 
this  mechanism  (formation  plus  migration  energy  of  the  pair)  was  computed  to 
be  3.9  eV,  in  excellent  agreement  with  experimental  data. 

iii)  In  the  heavily  doped  limit,  Mathiot  and  Pfister^  proposed  that  the  same 
mechanism  operates,  but  the  activation  energy  is  lowered  by  the  presence  of 
nearby  As  atoms.  In  particular,  during  migration,  as  the  vacancy  loops  around  a 
six-member  ring,  it  is  likely  to  pass  by  another  As  atom.  Mathiot  and  Fister 
recognized,  however,  that  such  proximity  with  another  As  atom  would  require 
extremely  high  concentrations.  They  resolved  this  quandary  by  invoking 
percolation  theory,  according  to  which  an  infinite  percolation  network  with  very 
high  As  concentrations  forms  when  the  total  concentration  exceeds  ~3x10^°  cm" 

The  activation  energy  for  the  mechanism  proposed  by  Mathiot  and  Pfister  was 
computed  to  be  ~1  eV  lower  than  the  simple  AsV  pair  mechanism.  It  would 
appear,  therefore,  that  the  percolation  model  can  account  for  the  mobility  of  As 
atoms  leading  to  clustering  at  moderate  temperatures  and  the  observed 
enhanced  diffusion  at  high  concentrations.  The  model,  however,  does  have  a 
problem  in  that  enhanced  diffusion  within  the  percolation  network  would  quickly 
lead  to  the  formation  of  AsnVm  clusters  that  would  break  up  the  continuity  of  the 
network.  In  contrast,  the  experiments  show  that  the  clustering  process  persists 
for  long  times. 

iv)  Calculations  reported  demonstrated  that  AS2V  complexes  are  mobile  and  can, 
therefore,  mediate  As  diffusion.  The  net  activation  energy  was  computed  to  be 
2.8  eV,  i.e.  essentially  the  same  as  in  the  mechanism  advocated  by  Mathiot  and 
Pfister.  The  key  difference  is  that  diffusion  mediated  by  AsaV  complexes  does 


not  need  a  continuous  percolation  network  but  only  high-concentration  regions. 
Such  regions  may  exist,  for  example,  when  the  average  concentration  is  just 
below  the  percolation  threshold  (pre-percolation  patches)  or  after  an  initial  rapid 
diffusion  for  average  concentrations  that  exceed  the  percolation  threshold  (post¬ 
percolation  patches).  In  the  latter  case,  the  initial  diffusion  leads  to  clustering 
that  breaks  up  the  percolation  network,  but  AsaV  complexes  continue  to  mediate 
diffusion  by  traveling  both  through  and  between  patches. 

V)  Mobile  AsaV  complexes  account  for  all  the  observed  phenomena  in  a  systematic 
way  (see  Appendix  A  for  more  details).  In  addition,  they  lead  to  predictions  on 
the  nature  of  the  dominant  complexes  in  deactivated  samples.  The  complexes 
AsV,  AS2V  and  AssV  are  the  only  one  which  can  break  up  into  an  isolated 
substitutional  As  and  a  mobile  complex.  In  each  case,  the  energy  to  “peel  off'  an 
As  atom  is  -1.5  eV,  which  is  the  same  as  the  value  of  the  reactivation  activation 
energy  Schwenker  et  al.^  extracted  from  their  experimental  data.^  Since  EXAFS 
data  established  that  in  deactivated  samples  virtually  all  As  atoms  are  second 
neighbors  to  each  others,  we  conclude  that  the  dominant  complexes  in 
deactivated  samples  are  AS2V  and  AS3V.  Higher-order  complexes  also  exist  and 
are  the  ones  that  persist  during  reactivation  at  1000°C  because  the  energy 
needed  to  break  them  up  is  much  larger  than  1 .5  eV. 

3.  Segregation  of  As  atoms  in  Si  grain  boundaries 

This  is  worked  carried  out  by  Pantelides  in  collaboration  with  A.  Maiti,  S.  J. 
Pennycook  and  M.  Chisholm  at  Oak  Ridge  National  Laboratory  and  is  described  in 
detail  in  Appendix  B,  where  references  to  original  literature  can  be  found.  The  issue 
is  that  virtually  99%  of  As  in  polycrystalline  Si  segregates  in  grain  boundaries.  This 
is  undesirable  because  As  in  grain  boundaries  is  electrically  inactive  and  one 
normally  would  like  to  dope  polycrystalline  Si  very  heavily  to  be  used  as  contact. 

The  segregation  energy  (defined  to  be  the  energy  difference  per  As  atom 
between  As  atoms  in  the  grain  boundary  and  the  bulk  crystal)  is  measured  to  be 
-0.5  eV  and  is  responsible  for  the  observed  behavior.  It  is  worthwhile  to  determine 
if  this  large  segregation  energy  is  an  intrinsic  property  of  grain  boundaries  in  Si  or  it 
is  a  consequence  of  defects  or  steps  where  As  atoms  bind  preferentially.  This 
question  became  particularly  relevant  when  it  was  determined  by  both  theory  and 
experiment  that  tilt  grain  boundaries  in  Si  without  steps  or  impurities  relax  to 
configurations  where  all  Si  atoms  are  fourfold-coordinated. 

Recent  theoretical  work  by  Arias  and  Joannopoulos^  found  that  As  atoms 
placed  at  fourfold-coordinated  sites  in  a  fully  relaxed  S=5  grain  boundary  in  Ge 
remain  fourfold  and  lower  their  energy  by  a  minimal  amount  over  bulk  sites  (-0.1 
eV).  They  suggested  that  the  observed  large  segregation  energy  is  likely  to  be 
caused  by  defects  or  steps. 

The  main  objective  of  our  work  was  to  explore  the  possibility  that  As  atoms  can 
be  incorporated  in  threefold  coordination  in  othenwise  perfect  grain  boundaries.  It 
was  found  that  threefold  coordination  is  possible  if  As  atoms  are  incorporated  in 
pairs.  The  phenomenon  is  quite  intriguing.  Two  nearest-neighbor  substitutional  As 


atoms  repel  each  other  trying  to  achieve  threefold  coordination  (their  preferred 
chemical  valence)  at  a  cost  of  elastic  energy  in  the  distorted  back  bonds.  The 
calculations  show  that  in  the  bulk  crystal  and  energy  gain  and  cost  roughly  balance 
out.  In  the  grain  boundary,  however,  the  cost  is  a  bit  less.  The  net  result  is  that  a 
dimer  binds  in  the  grain  boundary  through  repulsion!  The  segregation  energy  of 
such  dimers  was  found  to  be  in  the  range  0.1 -0.2  eV.  For  chains  of  dimers, 
however,  the  segregation  energy  increased  to  0.35  eV  (for  a  case  where  one 
actually  perfectly  ordered  As  atoms  with  equal  As-As  distances)  and  to  0.5  eV. 

The  dimer  calculations  were  done  for  two  distinct  S=5  boundary  structures,  both 
of  which  have  all  Si  atoms  fourfold  coordinated,  with  total  energies  that  are  very 
close  (one  is  higher  by  0.15  eV/atom).  The  As  chains  with  the  highest  segregation 
energy  are  in  fact  in  the  grain  boundary  structure  with  the  slightly  higher  energy. 
This  result  leads  to  the  intriguing  suggestion  that  As  atoms  segregating  in  a  grain 
boundary  may  actually  induce  a  structural  transformation  of  the  entire  grain 
boundary. 

In  addition  to  theoretical  investigations  of  As  impurities  in  Si  grain  boundaries, 
an  experimental  search  was  undertaken  by  Chisholm  and  Pennycook  for  evidence 
of  As  segregation.  Arsenic  was  implanted  and  annealed  for  a  prolonged  period  of 
time.  Micrographs  obtained  with  high-resolution  transmission  microscopy  using  the 
Z-contrast  technique  reveal  evidence  of  chains  of  As  atoms. 

Finally,  calculations  have  been  performed  to  investigate  whether  the  ordering  of 
As  impurities  in  grain  boundaries  occurs  also  in  isolated  dislocations.  The  results 
indicate  that  chains  of  As  dimers  again  form  along  the  dislocation  pipe  with  very 
large  segregation  energies  (~1  eV).  These  results  are  described  in  detail  in 
Appendix  C. 

4.  Oxygen  migration  in  Si 

This  is  a  case  that  on  the  face  appears  very  simple,  yet  turns  out  to  hide 
enormous  complexity.  The  work  was  carried  out  by  Ramamoorthy  and  Pantelides 
and  is  described  in  detail  in  Appendix  D,  where  references  to  the  original  literature 
can  be  found. 

Oxygen  has  long  been  known  to  occupy  a  bridge  site  between  neighboring  Si 
atoms,  forming  a  buckled  Si-O-Si  chain  virtually  identical  to  the  Si-O-Si  chains  in 
Si02.  Its  migration  has  been  measured  over  a  range  of  more  than  1000  degrees 
and  has  been  found  to  have  a  simple  Arrhenius  behavior  with  a  2.5-eV  activation 
energy.  The  migration  ought  to  be  a  simple  process:  the  0  atom  jumps  from  one 
bridge  site  to  the  next.  The  saddle  point  should  be  at  the  midpoint  of  the  O  path. 

Theory  over  the  years  has  stumbled  on  this  one.  Calculations  of  the  activation 
energy  have  ranged  from  1.2  eV  to  4.1  eV.  Some  authors  found  2.5  eV.  Then,  a 
paper  by  Jiang  and  Brown  in  1995  reported  calculations  using  classical  interatomic 
potentials.  These  authors  stepped  the  O  atom  along  its  entire  path,  relaxing  the  Si 
atoms  at  each  step.  They  found  that  the  maximum  energy,  i.e.  the  saddle  point  was 
at  a  point  after  the  O  atom  passed  the  midpoint  of  its  path.  This  saddle  point  energy 


was  2.5  eV.  However,  Jiang  and  Brown"*  did  not  notice  that  the  total-energy  profile 
they  computed  along  the  path  did  not  satisfy  a  basic  symmetry  requirement,  namely 
time  reversal,  which  requires  that  the  energy  profile  be  symmetric  about  the 
midpoint  of  the  path. 

We  performed  systematic  first-principles  calculations  and  unveiled  an  enormous 
amount  of  complexity  while  accounting  for  both  the  experimental  data  and  the 
earlier  theoretical  work.  It  was  found  that  the  motion  of  the  O  atom  is  intimately 
coupled  to  the  motion  of  the  “central  Si  atom”  -  the  one  that  participates  in  both  the 
initial  and  the  final  Si-O-Si  chain.  Both  these  atoms  face  barriers.  In  order  to 
describe  the  migration  process  one  needs  a  total-energy  hypersurface  in  several 
dimensions  that  reveals  a  multiplicity  of  paths  and  a  “saddle  ridge”  with  a  height  of 
~  2.5  eV.  Details  and  figures  can  be  found  in  Appendix  D. 

More  recent  work  (Appendix  E)  examined  the  effect  of  H  on  O  diffusion.  It  was 
found  that  H  also  occupies  a  buckled  bridge  position  next  to  the  O.  Migration  of  the 
O  atom  now  proceeds  as  follows:  the  O  atom  heads  toward  the  bridge  site 
occupied  by  H;  the  latter  moves  toward  the  next  bridge  site  in  front  of  it.  This 
synergism  led  to  a  lowering  of  the  activation  by  about  1  eV,  in  agreement  with 
observations. 

5.  Atomic-scale  simulation  of  heat  flow 


Most  molecular  dynamics  simulations  are  carried  out  within  one  of  the 
equilibrium  ensembles  of  statistical  mechanics.  In  order  to  carry  out  simulations  of 
a  nonequilibrium  process  such  as  heat  flow  one  has  to  meet  significantly  more 
stringent  requirements.  We  have  pursued  heat  flow  calculations  by  maintaining  a 
difference  in  temperature  at  the  two  ends  of  a  slab  of  Si  and  demonstrated  that  by 
using  long  simulations  it  is  possible  to  establish  local  equilibrium  in  thing  layers  and 
set  up  a  smooth  temperature  gradient.  We  pursued  such  calculation  for  a  grain 
boundary  in  Si  using  classical  interatomic  potentials  and  demonstrated  that  it  is 
possible  to  compute  the  corresponding  Kapitza  resistance.  Details  are  given  in 
Appendix  F. 

6.  Mesoscopic  calculations  of  void  growths 

Electrical  current  is  known  to  induce  atomic  diffusion  which  leads  to  void  growth 
or  extrusions  that  can  sever  or  short  metal  interconnects  in  integrated  circuits.  Our 
work  focused  on  void  growth.  Results  are  described  in  Appendix  G. 
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Several  complex  dynamical  phenomena  have  been  observed  in  heavily  doped  Si,  but  a  comprehensive 
account  of  the  underlying  atomic-scale  processes  is  lacking.  We  report  a  wide  array  of  first-principles 
calculations  in  terms  of  which  we  give  such  a  comprehensive  account.  In  particular,  we  find  that 
vacancies  (V),  AsV  pairs,  AS2V'  complexes,  and  higher-order  As„V„  complexes  play  distinct  roles  in 
the  observed  dopant  deactivation,  reactivation,  and  anomalous  diffusion.  The  latter  is  mediated  by 
mobile  As2  V  complexes  that  form  in  “prepercolation”  patches  of  a  very  high  dopant  concentrations  and 
gives  rise  to  fast  As  clustering  at  moderate  temperatures.  Our  results  are  quantitative  and  in  agreement 
with  experimental  numbers  where  available.  [80031-9007(96)00428-0] 
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For  some  time  now  it  has  been  possible  to  prepare  Si 
specimens  with  very  high  concentrations  of  dopants  such 
as  As  (as  high  as  10''  cm^;  i.e.,  roughly  one  As  atom  per 
50  host  atoms)  with  a  corresponding  high  conductivity. 
Such  specimens  exhibit  a  wide  array  of  complex  dynamical 
phenomena  that  have  been  investigated  extensively  over 
the  last  25  years.  These  investigations  fall  into  two  distinct 
categories. 

(i)  Beginning  in  the  1 970s,  a  number  of  investigators  [1  — 
4]  found  that  annealing  of  heavily  doped  samples  at  mod¬ 
erate  temperatures  (400-500  °C)  results  in  a  dramatic  drop 
in  conductivity.  The  phenomenon  has  been  universally  at¬ 
tributed  to  the  formation  of  electrically  inactive  As  com¬ 
plexes.  This  deactivation  of  dopants  is  partially  reversible 
by  heating  to  higher  temperatures  (800- 1 000  °C).  The  ma¬ 
jor  focus  of  most  investigations  has  been  the  identification 
of  the  dominant  electrically  inactive  complex.  There  is  ex¬ 
tensive  experimental  evidence  that  the  dominant  defect  is 
an  As„  complex,  where  V  stands  for  a  vacancy,  but  there 
is  no  consensus  on  the  value  of  n  [5-11].  TTie  most  per¬ 
plexing  feature  of  the  observations  is  that  at  400-500  ®C, 
in  lightly  doped  Si,  As  diffusion  is  negligible.  Yet,  some¬ 
how,  at  these  temperatures,  in  heavily  doped  specimens, 
the  As  atoms  are  rapidly  rearranged  [1-3]  from  isolated 
substitutional  sites  to  form  As„  V  clusters  [9-11]  and  pre¬ 
cipitates  [12],  generating  silicon  self-interstitials  [11].  No 
atomistic  mechanism  for  the  deactivation  and  subsequent 
reactivation  has  been  proposed  so  far. 

(ii)  A  different  set  of  investigators  studied  As  diffu¬ 
sion  in  heavily  doped  samples  at  temperatures  in  the  range 
of  1000-1 100  °C.  Early  investigations  [5]  showed  that 
the  As  diffusion  constant  increased  linearly  with  dopant 
concentration  «  for  n  <  3  X  10“  cm"^.  For  n  >  3  X 
10^®  cm  ^  the  diffusion  constant  decreased.  This  was  as¬ 
cribed  to  clustering  [5].  However,  later  rapid  thermal  an¬ 
nealing  (RTA)  experiments  [13,14]  revealed  increased  As 
diffusivity  for  n  >  3  X  10*°  cm  “I  Other  RTA  studies 
[15]  of  the  diffusion  of  small  concentrations  of  As  in  heav¬ 
ily  P-doped  Si  also  showed  a  dramatic  increase  in  As  dif- 

0031-9007/96/76(25)/4753(4)$10.00 


fusivity  for  donor  concentrations  above  ~2  X  10^®cm~^. 
The  most  intriguing  suggestion  to  account  for  the  observed 
enhanced  dopant  diffusion  above  a  critical  doping  thresh¬ 
old  is  the  formation  of  a  percolation  network  of  dopant 
atoms  [15,16].  However,  other  theories  [17]  account  for 
the  enhanced  diffusion  without  invoking  percolation.  No 
attempt  has  been  made  thus  far  to  give  a  unifying  interpre¬ 
tation  of  the  diffusion  data  and  the  deactivation  or  reacti¬ 
vation  data. 

In  this  Letter,  we  report  a  comprehensive  set  of  first- 
principles  calculations  in  terms  of  which  we  obtain  a 
detailed  description  of  the  dynamical  atomic-scale  pro¬ 
cesses  that  occur  in  heavily  As-doped  Si.  These  re¬ 
sults  allow  us  to  give  a  systematic  account  of  both  the 
deactivation-reactivation  observations  and  the  anomalous 
diffusion  data.  In  particular,  we  find  that  vacancies,  AsV 
pairs,  AsiV  complexes,  and  higher-order  As„V,„  com¬ 
plexes  play  distinct  roles  in  the  observed  deactivation,  re¬ 
activation,  and  anomalous  diffusion  phenomena.  A  per¬ 
colation  network,  which  in  principle  may  exist  in  heavily 
doped  specimens  prior  to  any  thermal  treatments,  would 
in  fact  quickly  break  up  due  to  the  formation  of  As„Vm 
complexes  and  thus  would  not  sustain  long-term  As  dif¬ 
fusion  and  clustering.  We  find  that  fast  As  diffusion  at 
high  doping  levels  is  mediated  by  mobile  AS2V complexes 
that  form  in  “prepercolation”  patches  of  very  high  dopant 
concentrations  and  give  rise  to  rapid  As  clustering  at  mod¬ 
erate  temperatures.  Our  predictions  are  quantitative  and  in 
agreement  with  experimental  numbers,  where  available. 

The  calculations  were  based  on  the  density  functional 
method,  using  the  local  density  approximation.  Pseudopo¬ 
tentials  [18]  were  used  for  Si  and  As.  Electronic  wave 
functions  were  expanded  in  a  plane  wave  basis,  with  a 
15  Ry  cutoff.  Tests  for  bulk  Si  gave  a  lattice  parame¬ 
ter  of  5.42  A  and  a  bulk  modulus  of  0.99  Mbar.  Defect 
configurations  were  simulated  using  the  supercell  method. 
The  formation  energies  of  individual  defects  changed  by 
only  0.25  eV  on  going  from  16  to  32  site  supercells,  in¬ 
dicating  that  the  results  from  32  site  supercells  were  well 
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converged.  Two  special  k  points  were  used  [19].  Migra¬ 
tion  energies  of  the  As- vacancy  complexes  were  calculated 
using  64  site  supercells.  Calculated  forces  on  atoms  were 
used  to  relax  each  structure  to  equilibrium. 

First,  we  report  our  calculations  of  the  formation  en¬ 
ergies  of  defect  clusters  containing  As  atoms  and  native 
point  defects  in  Table  I.  Our  results  show  that  As  atoms 
have  negligible  affinity  for  each  other  but  bind  strongly 
to  vacancies  and  interstitials.  The  binding  energy  per 
As  atom  in  a  cluster  containing  a  vacancy  is  apprecia¬ 
bly  larger  than  that  containing  an  interstitial,  reflecting 
the  preference  of  As  for  threefold  coordination.  In  fact, 
each  time  an  As  atom  is  taken  from  an  isolated  substi¬ 
tutional  fourfold-coordinated  site  and  placed  next  to  a  va¬ 
cancy,  there  is  an  energy  gain  of  about  1.5  eV.  All  As„ 
complexes  with  n  >  2m  have  negative  formation  ener¬ 
gies,  including  the  AS4V  complex,  where  our  result  is 
in  agreement  with  the  calculation  reported  by  Pandey  et 
ai  [9].  These  authors  studied  only  the  As4y  complex 
and  asserted  that,  because  of  its  negative  formation  energy, 
it  is  the  dominant  defect  in  deactivated  samples.  How¬ 
ever,  as  we  see,  this  complex  is  only  one  among  a  hier¬ 
archy  of  such  complexes.  In  these  results,  we  observe 
a  clear  trend:  At  any  doping  level  compared  with  iso¬ 
lated  fourfold-coordinated  configurations,  thermodynam¬ 
ics  would  favor  the  agglomeration  of  As  atoms  with  va¬ 
cancies  in  the  form  of  As complexes  or  the  attach¬ 
ment  of  As  atoms  to  microvoid  surfaces.  In  other  words. 
As-doped  Si  is  metastable  at  any  doping  level  and  the 
preparation  of  electrically  active  samples  is  possible  only 
because  of  kinetic  limitations. 

The  results  presented  above  allow  us  to  conclude  that 
the  defects  that  underlie  deactivation  and  reactivation  of 
heavily  doped  Si  are  As„  complexes.  This  conclusion 
is  fully  consistent  with  earlier  analysis  of  experimental  data 
[6,8-11].  It  is  clear,  however,  that  formation  energies 


TABLE  I.  The  formation  energies,  in  eV,  of  the  complexes  of 


As  atoms  with  a  vacancy. 
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Formation 

energy 

Binding  energy 
per  As  atom 
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3.78 
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1.31 

AS2V 

0,82 
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-0.54 

1.65 

AS6V2 

-3.23 

1.57 

7  (interstitial) 

4.14 

AS2/ 

3.76 
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FIG.  1.  Migration  of  AsV.  Si:  white  circle,  As:  dark  circle, 
V:  shaded  circle,  (a)  Initial  configuration;  (b)  As  and  V 
have  exchanged  positions;  (c)  saddle  configuration:  V  is  third 
neighbor  of  As;  (d)  final  configuration;  one  diffusion  step  is 
complete. 


alone  (thermodynamics)  do  not  tell  us  what  role  each 
As„  Vm  complex  plays,  which  one  is  dominant  if  any,  what 
are  the  mechanisms  that  allow  them  to  form  only  in  heavily 
doped  materials,  and  how  diffusion  is  enhanced.  We  must, 
therefore,  investigate  the  dynamical  processes  (kinetics) 
that  are  possible  only  in  heavily  doped  materials. 

It  is  helpful  to  start  with  the  process  that  underlies  As 
diffusion  in  lightly  doped  Si.  It  is  well  known  that  dif¬ 
fusion  proceeds  primarily  by  the  formation  and  migration 
of  AsV  pairs  [20].  The  migration  of  AsV  is  illustrated  in 
Fig.  1:  The  vacancy  loops  around  the  As  atom  for  the  pair 
to  make  one  step.  The  diffusion  activation  energy  of  the 
AsV  complex  is  the  sum  of  its  formation  energy  and  mi¬ 
gration  energy.  From  Table  I,  the  formation  energy  of  the 
AsV  complex  is  2.5  eV.  The  migration  energy  is  made  up 
of  two  contributions:  the  change  in  energy  when  the  va¬ 
cancy  goes  from  the  nearest  neighbor  to  the  third  neighbor 
position  of  the  As  atom  and  the  energy  of  vacancy  mi¬ 
gration,  which  we  calculate  to  be  0.9  and  0.5  eV,  respec¬ 
tively.  This  gives  a  net  activation  energy  for  AsV  diffusion 
of  3.9  eV,  in  good  agreement  with  experiment  [20].  With 
such  a  large  activation  energy,  negligible  As  migration  oc¬ 
curs  at  moderate  temperatures. 

In  the  heavily  doped  limit,  Mathiot  and  Pfister  (MP)  [16] 
have  suggested  that  enhanced  diffusion  comes  about  as  fol¬ 
lows:  AsV  pairs  during  their  migration  are  not  isolated. 
Specifically,  MP  considered  the  case  there  two  As  atoms 
are  fifth  neighbors.  If  one  of  these  As  atoms  is  connected 
to  a  vacancy,  then,  during  the  migration  of  the  AsV  pair, 
when  the  V  is  a  third  neighbor  of  one  As  atom,  it  would 
be  a  second  neighbor  of  the  other.  The  net  effect  is  a  low¬ 
ering  of  the  migration  energy.  For  a  uniform  As  distribu¬ 
tion,  such  a  cluster  would  not  occur  even  at  an  As  con¬ 
centration  as  high  as  10^^  cm”^,  for  As  atoms  would  be 
eighth  neighbors,  on  average.  MP  proposed  a  resolution  to 
this  quandary  by  invoking  percolation  theory.  They  found 
that  when  the  As  concentration  exceeds  *^3  X  10^^  cm”^, 
fluctuations  in  the  As  distribution  would  result  in  the  pres¬ 
ence  of  an  infinite  network  of  As  atoms  in  the  sample, 
in  which  any  two  As  atoms  are  fifth  neighbors  of  each 
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other  (or  closer).  They  described  a  scenario  for  enhanced 
As  diffusion  above  this  critical  doping  threshold,  due  to 
accelerated  diffusion  of  vacancies  through  this  “infinite 
percolation  network.*’  We  have  performed  the  relevant 
calculations  in  this  limit  using  supercells  and  found  that 
both  the  formation  and  migration  energy  of  an  AsV  pair 
are  lowered,  leading  to  a  decrease  in  the  diffusion  acti¬ 
vation  energy  of  about  1  eV.  This  is  consistent  with  the 
number  obtained  by  MP  [16]  from  simulations.  A  sim¬ 
ple  calculation,  assuming  no  significant  change  in  the  pre¬ 
exponential  factor  in  the  expression  for  the  diffusion  co¬ 
efficient  [D  =  DoCxp(-Q/kbT)],  indicates  that  the  rate 
of  As  diffusion  in  a  heavily  doped  specimen  at  '^450  °C 
is  comparable  to  that  in  a  lightly  doped  specimen  at 
~900  °C.  This  can  explain  how  fast  As  migration  can  oc¬ 
cur  at  moderate  temperatures. 

However,  there  is  a  problem  with  the  above  scenario. 
The  percolation  model  yields  a  quantitative  account  of 
enhanced  diffusion  for  As  concentrations  above  -^3  X 
10^^  cm”^,  but  it  is  not  consistent  which  the  fact  that 
enhanced  diffusion  causes  As  clustering  with  vacancies. 
Such  clustering  would  break  up  the  percolation  network 
into  high  concentration  patches  and  halt  enhanced  diffu¬ 
sion  and  further  clustering  after  a  few  seconds  of  annealing. 
This  would  be  in  contradiction  with  deactivation  observa¬ 
tions  [1,2]  that  indicate  continued  As  clustering  over  many 
hours  at  moderate  temperatures. 

We  note  that  other  workers  have  claimed  to  account  for 
enhanced  diffusion  [17]  in  heavily  doped  specimens  with¬ 
out  invoking  percolation.  Here,  we  propose  an  atomistic 
mechanism  that  does  not  require  a  fully  connected  perco¬ 
lation  cluster  and  yet  explains  the  diffusion  data  [5,14], 
namely,  a  mobile  As2y  complex.  The  formation  energy  of 
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FIG.  2.  Migration  of  As  2V  (a)  Initial  configuration;  (b)  one 
As  and  V  have  exchanged  positions;  (c),(d)  translation;  (e),(f) 
reorientation  starting  with  configuration  (b). 


AS2V,  from  Table  I,  is  0.8  eV.  Its  migration  is  illustrated 
schematically  in  Fig.  2.  Two  basic  steps  are  depicted:  the 
sequence  (abed)  and  the  sequence  (abef ),  with  comparable 
migration  energies.  The  net  diffusion  activation  energy  is 
2.7  eV  for  the  As2  V  complex,  which  is  1.2  eV  lower  than 
that  of  the  As  V  pair.  This  is  virtually  identical  to  the  result 
we  obtained  for  AsV  migration  in  the  percolation  network. 
Thus,  based  on  activation  energies  alone,  one  cannot  dis¬ 
tinguish  between  the  percolation  model  and  the  migration 
of  As 2  V.  While  earlier  investigators  [5,7,10]  had  proposed 
that  AS2V'  complexes  [21]  are  present  in  deactivated  spec¬ 
imens,  they  were  not  considered  mobile.  In  modeling  dif¬ 
fusion  of  Sb  in  P-doped  Si,  Larsen  and  Weyer  [22]  had 
considered  diffusion  mediated  by  SbPV  complexes,  but  in 
later  work  [15]  abandoned  this  mechanism  in  favor  of  the 
percolation  model.  We  believe,  however,  that  As 2V plays 
a  key  role  in  enhanced  diffusion  in  heavily  doped  Si  as 
it  does  not  need  a  continuous  percolation  network,  some¬ 
thing  that  would  not  be  sustained  when  As  agglomera¬ 
tion  into  As„Vm  complexes  occurs.  We  propose  that  en¬ 
hanced  As  diffusion  starts  below  the  percolation  thresh¬ 
old,  calculated  by  MP  [16]  to  be  3  X  10^^  cm”^.  High- 
concentration  “prepercolation  patches”  in  which  two  As 
atoms  are  fifth  neighbors  or  closer  enable  the  formation 
of  large  numbers  of  AS2V  complexes.  These  AS2V  com¬ 
plexes  can  migrate  fast  through  the  sample,  even  at  mod¬ 
erate  temperatures.  During  their  motion,  they  react  with 
other  defects  to  form  larger,  immobile  complexes.  Model¬ 
ing  diffusion  data  on  the  basis  of  these  reactions  is  beyond 
the  scope  of  this  Letter. 

We  now  summarize  the  main  experimental  observa¬ 
tions  and  how  the  theory  described  above  accounts  for 
them. 

(1)  Diffusion  data. — Early  observations  [5]  showed  that 
the  diffusion  constant  of  As  decreased  for  dopant  concen¬ 
trations  greater  than  3  X  10^^  cm”^.  However,  later  RTA 
experiments  [14]  on  As-implanted  Si  found  fast  As  redis¬ 
tribution  (within  20  s)  which  reduced  the  peak  As  concen¬ 
tration  to  (2.5-3)  X  10^*^  cm”^,  if  the  as-implanted  peak 
As  concentration  exceeded  this  value.  Otherwise,  negli¬ 
gible  As  migration  was  observed.  RTA  experiments  [15] 
on  small  concentrations  of  As  in  heavily  P-doped  Si  found 
enhanced  diffusion  for  donor  concentrations  greater  than 
—2.0  X  10^^  cm  The  early  diffusion  studies  were  car¬ 
ried  out  over  several  hours.  The  decreased  As  diffusivity 
at  high  doping  levels  was  correctly  attributed  to  As  clus¬ 
tering  [5].  The  fast  As  migration  observed  in  RTA  experi¬ 
ments  is  due  to  As 2V  migration  over  a  very  short  time  scale 
(6-60  s),  during  which  extensive  clustering  does  not  oc¬ 
cur.  The  investigations  of  As  diffusion  in  heavily  P-doped 
Si  establish  that  the  onset  of  enhanced  dopant  diffusion  is 
indeed  below  the  percolation  threshold.  The  enhanced  As 
diffusion  observe  din  these  latter  experiments  is  mediated 
by  AsPV  complexes.  The  extracted  activation  energy  [15] 
of  2.7  eV  is  identical  to  that  of  the  As  2V  complex. 

(2)  Deactivation  data. — During  deactivation  at  moder¬ 
ate  temperatures  (400-600  ®C),  the  rate  of  deactivation  is 
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highest  initially  and  then  continues  at  a  slower  rate  [1,2,4]. 
Initially  As  atoms  are  close  to  each  other  in  prepercolation 
patches.  Migrating  As2  V  complexes  have  a  high  cross  sec¬ 
tion  for  reaction  inside  such  patches,  giving  rise  to  larger 
immobile  complexes.  Thus,  the  initial  deactivation  rate  is 
high.  Once  a  degree  of  clustering  takes  place,  As  atoms  are 
further  from  each  other,  on  average,  and  As  2  V  complexes 
have  to  travel  greater  distances  to  cause  further  deactiva¬ 
tion.  At  some  point,  deactivation  becomes  extremely  slow 
and  effectively  stops.  The  distribution  of  complexes  is  de¬ 
termined  by  kinetics,  for  the  system  does  not  reach  ther¬ 
modynamic  equilibrium  at  these  temperatures. 

(3)  Reactivation  data. — Annealing  of  deactivated  spec¬ 
imens  at  elevated  temperatures  (800-1000  °C)  is  found  to 
reactivate  carriers.  Schwenker,  Pan,  and  Lever  [1],  as¬ 
suming  that  one  complex  is  predominant  in  deactivated 
specimens,  obtained  the  activation  energy  for  the  reacti¬ 
vation  to  be  1 .6  eV  and  the  number  of  As  atoms  in  the 
complex  to  be  between  1.7  and  2.8  As  atoms.  For  con¬ 
centrations  of  order  2  X  10^®  cm“^,  most  of  the  initial  As 
concentration  is  reactivated.  Reactivation  is  a  result  of  the 
breakdown  of  complexes  formed  at  lower  temperatures, 
releasing  As  atoms,  and  requires  the  removal  of  an  AsV 
pair  or  an  As 2  V  complex,  for  these  are  the  mobile  species. 
From  Table  I,  we  find  that  the  energy  required  to  break  up 
AsV  (AsV  As  -h  As),  AS2V  (AS2Y  AsV  +  As), 
and  AS3V  (AS3V  — ►  AS2V  +  As),  are  all  about  1.5  eV, 
in  excellent  agreement  with  the  data.  Larger  complexes 
do  not  decompose  into  As  +  Y,  as  Y  is  not  mobile.  They 
would  have  to  break  up  into  larger  fragments  and  that  re¬ 
quires  more  energy.  Therefore,  we  conclude  that  AsV, 
AS2K  and  AS3  V  may  be  present  in  deactivated  specimens. 
AsV  is  ruled  out  by  experiments  [9,11]  which  find  that  most 
As  atoms  in  deactivated  specimens  have  another  As  atom 
as  a  second  neighbor.  The  fact  that  most  of  the  As  atoms 
were  reactivated  at  a  doping  level  of  2  X  10^^  cm”^  indi¬ 
cates  that  As  2V  and  As  3  V  are  the  dominant  complexes  in 
deactivated  specimens,  near  the  enhanced-diffusion  thresh¬ 
old. 

(4)  The  saturation  of  the  carrier  concentration. — The 
carrier  concentration  in  heavily  doped  Si  does  not  exceed 

X  10^®  cm“^  at  elevated  temperatures  (1000 ®C)  [3], 
even  if  the  total  As  concentration  exceeds  this  value.  This 
concentration,  as  MP  [16]  showed,  marks  the  threshold 
above  which  an  infinite  percolation  network  of  As  atoms 
is  present  initially  in  the  sample.  If  the  initial  As  concen¬ 
tration  exceeded  this  value,  then,  in  the  very  early  stages 
of  annealing,  a  very  high  flux  of  vacancies  would  flood 
the  sample  along  the  percolation  network.  Thus,  the  As 
concentration  in  excess  of  the  critical  value  would  become 
bound  very  rapidly  in  As„V;„  complexes  with  n  >  2m 
(like  AS4V).  Such  complexes  would  not  decompose  even 
at  1000  "C. 

In  conclusion,  we  reported  extensive  calculations  of 
formation  energies,  migration  energies,  and  defect  reaction 


energies  in  terms  of  which  we  have  been  able  to  give 
a  systematic  account  of  the  major  observations  of  both 
deactivation  or  reactivation  and  enhanced  diffusion  in 
heavily  doped  Si. 
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Recent  theoretical  work  found  that  isolated  As  impurities  in  Ge  grain  boundaries  exhibit  minimal 
binding,  leading  to  the  suggestion  that  the  observed  segregation  is  likely  to  occur  at  defects  and 
steps.  We  report  ab  initio  calculations  for  As  in  Si  and  show  that  segregation  is  possible  at  defect- 
free  boundaries  through  the  cooperative  incorporation  of  As  in  threefold-coordinated  configurations: 
As  dimers,  or  ordered  chains  of  either  As  atoms  or  As  dimers  along  the  grain  boundary  dislocation 
cores.  Finally,  we  find  that  As  segregation  may  drive  structural  transformations  of  certain  grain 
boundaries.  [S003 1  -9007(96)008 1 1  -3] 


PACS  numbers:  61.72.Mm,  68.35.Dv,  68.55.Ln,  71.15.Pd 

Polycrystalline  semiconductors  are  used  in  microelec¬ 
tronics.  Dopants,  which  control  the  electrical  properties 
of  these  materials,  are  known  to  segregate  in  grain  bound¬ 
aries  in  electrically  inactive  configurations  [1-4].  The 
fraction  of  dopants  in  the  grain  boundaries  is  governed  by 
the  segregation  energy,  i.e.,  the  energy  difference  between 
a  dopant  atom  in  the  grain  boundary  and  a  dopant  atom  in 
the  bulk  crystal.  For  the  specific  case  of  As  segregation 
in  Si  and  Ge  grain  boundaries,  experimental  values  of  the 
segregation  energies  range  from  0.41  to  0,65  eV  [1-4]. 

Clearly  it  would  be  very  important  for  technologists  to 
know  if  the  large  segregation  energy  is  an  intrinsic  property 
of  a  defect-free  grain  boundary  or  is  caused  by  defects 
that  could,  in  principle,  be  avoided.  The  origin  of  the 
segregation  energy  has  not,  however,  been  accounted  for 
so  far.  Experiment  [5]  and  theory  [6,7]  have  established 
that  tilt  grain  boundaries  in  undoped  Si  and  Ge  rebond 
so  that  all  host  atoms  are  fourfold  coordinated.  The 
only  first-principles  theoretical  study  of  segregation  was 
reported  recently  by  Arias  and  Joannopoulos  [8].  These 
authors  examined  the  segregation  energies  of  isolated  As 
atoms  placed  at  different  substitutional  sites  in  a  Ge  grain 
boundary  and  found  values  only  of  order  0.1  eV.  They 
proposed  that  the  observed  segregation  energies  are  likely 
to  arise  from  As  atoms  bound  to  steps  or  other  defects.  No 
calculations  were  pursued  to  explore  such  possibilities. 

The  motivation  for  the  present  work  was  the  recogni¬ 
tion  that  the  observed  large  segregation  energies  may  occur 
in  defect-free  grain  boundaries  because  As  atoms  achieve 
their  preferred  threefold  coordination,  as  they  are  known 
to  do  in  amorphous  Si  and  Ge  [9].  Simple  bond  counting 
suggests  that  if  a  single  As  atom  were  to  achieve  three¬ 
fold  coordination  in  a  Si  or  Ge  grain  boundary,  at  least 
one  Si  (Ge)  atom  would  have  to  have  odd  coordination 
(3  or  5),  which  is  energetically  costly.  Clearly  threefold- 
coordinated  As  atoms  would  be  far  more  likely  if  they 
were  incorporated  in  a  grain  boundary  in  a  cooperative 
manner,  at  least  two  at  a  time.  The  simplest  possibility 
would  be  two  As  atoms  at  nearest-neighbor  sites.  We  have 
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performed  systematic  first-principles  calculations  for  such 
possibilities  in  Si,  which  is  technologically  more  impor¬ 
tant,  and  found  that  As  dimers  do  form  in  grain  bound¬ 
aries  with  segregation  energies  up  to  0.2  eV  per  As  atom. 
Dimer  binding  occurs  through  repulsion:  The  two  As 
atoms  repel  each  other  seeking  to  achieve  threefold  coordi¬ 
nation  (Fig.  1);  the  energy  gain  from  this  chemical  rebond¬ 
ing  is  larger  than  the  elastic  energy  cost  from  backbond 
distortions.  Furthermore,  we  find  that  larger  segregation 
energies  (0.3 -0.5  eV  per  As  atom)  are  achieved  through 
the  formation  of  chains  of  As  dimers  or  fully  ordered 
chains  of  threefold-coordinated  As  atoms.  Finally,  we  find 
that  the  formation  of  such  chains  may  induce  a  structural 
transformation  of  certain  grain  boundaries. 

Calculations  were  performed  for  a  S  =  5{310}(001) 
symmetric  tilt  boundary  in  Si.  This  grain  boundary  is 
parallel  to  the  {310}  plane  of  the  original  crystalline  lat¬ 
tice.  It  has  a  minimum  periodicity  of  one  conventional 
lattice  parameter  {a  =  5.431  A)  in  the  (001)  direction  and 
a  periodicity  of  a^J5|2  —  8.587  A  in  the  direction  per¬ 
pendicular  to  the  (001)  axis.  In  the  discussion  below, 
we  follow  the  convention  that  the  jc  axis  is  perpendicu¬ 
lar  to  the  grain  boundary  plane,  the  z  axis  is  parallel  to 
the  (001)  axis  of  the  original  crystalline  lattice,  and  the 
y  axis  runs  parallel  to  the  grain  boundary  plane  in  a  di¬ 
rection  perpendicular  to  the  z  axis.  We  used  periodic  su¬ 
percells  that  contain  two  oppositely  oriented  S  =  5  grain 
boundaries.  In  a  supercell  with  N  planes  of  atoms  along 


Si  Si 


FIG.  1.  Schematic  showing  how  the  member  atoms  of  an 
As  dimer  placed  substitutionally  in  neighboring  Si  sites  would 
move  away  from  each  other  into  threefold  coordination,  and 
thus  lower  energy. 
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the  X  direction,  the  grain  boundaries  are  separated  by 
N/2  planes.  The  actual  supercells  used  in  the  calcula¬ 
tions  will  be  discussed  later. 

The  calculations  were  based  on  density  functional  the¬ 
ory  [10]  with  local  exchange-correlation  energy  as  pa¬ 
rametrized  by  Perdew  and  Zunger  [11].  The  atomic  cores 
are  represented  by  nonlocal,  norm-conserving  pseudopo¬ 
tentials  of  the  Kerker  type  [12]  in  a  separable  Kleinmann- 
Bylander  form  [13],  and  defined  on  a  real-space  grid  [14]. 
The  calculations  were  performed  using  the  code  CETEP 
[15],  which  was  run  on  128  processors  of  the  Intel  Paragon 
XP/S  35  at  Oak  Ridge  National  Laboratory.  The  inte¬ 
gration  over  the  Brillouin  zone  was  performed  using  two 
special  k  points  chosen  according  to  the  Monkhorst-Pack 
scheme  [16].  The  electronic  wave  functions  were  ex¬ 
panded  in  a  plane  wave  basis  set  with  an  energy  cut¬ 
off  of  150  eV,  verified  to  yield  accurate  lattice  constant 
and  bulk  modulus  for  the  pure  crystal.  For  each  geome¬ 
try  the  electronic  wave  functions  were  first  relaxed  by 
the  conjugate  gradient  scheme  of  Payne  et  aL  [17]  until 
they  reached  a  local  minimum  (the  Bom-Oppenheimer 
surface).  The  ions  were  then  moved  according  to  the 
Hellman-Feynman  forces  until  the  largest  force  on  any  ion 
in  any  direction  was  less  than  0.08  eV/A.  Energy  changes 
due  to  changes  in  supercell  dimensions,  known  to  be  small 
[8],  were  neglected. 

Calculations  were  performed  with  V  =  30  and 
40  planes  of  atoms,  with  the  corresponding  supercells 
containing,  respectively,  60  and  80  atoms  for  the  mini¬ 
mum  periodicity  in  the  grain  boundary  plane.  Computed 
segregation  energies  changed  by  less  than  0.02  eV,  indi¬ 
cating  that  the  two  grain  boundaries  in  the  supercell  are 
adequately  isolated  from  each  other.  In  order  to  isolate 
the  As  dimers  from  each  other  in  the  grain  boundary 
plane,  calculations  were  performed  with  double  the  primi¬ 
tive  cell  and  N  =  30.  Most  calculations  were  performed 
using  only  the  primitive  cell  in  the  grain  boundary  plane 
and  N  =  40,  resulting  in  chains  of  interacting  dimers. 

The  calculations  for  the  undoped  2  =  5  grain  boundary 
yielded  a  fully  relaxed  structure,  which  we  label  GBl, 
that  is  the  same  as  described  in  Ref.  [8]  for  Ge.  We 
found,  however,  a  second  low-energy  structure,  which 
we  label  GB2,  with  a  total  energy  that  is  higher  by  only 
0.15  eV  per  periodic  segment  of  the  grain  boundary  plane. 
Figures  2  and  3  display  xy  projections  of  the  structures. 
The  two  structures  differ  in  the  relative  z  shift  of  the 
two  grains  forming  the  boundary,  by  a/S  =  0.68  A,  and 
also  in  the  nature  of  the  dislocation  cores  comprising  the 
boundaries.  The  cores  of  GBl,  formed  by  terminating 
planes  coming  from  the  same  grain,  are  of  the  pure  edge 
type  (b  =  f<100».  The  cores  of  GB2  are  formed  by 
terminating  planes  coming  from  different  grains,  have  the 
Burger’s  vector  at  45°  to  the  z  direction  (b  =  f  (101)), 
and  therefore  have  mixed  screw  and  edge  character. 

We  studied  segregation  of  As  atoms  in  both  GBl  and 
GB2  because  the  two  structures  contain  different  dislo¬ 
cation  cores  that  are  components  of  many  different  grain 


FIG.  2.  2D  projection  (normal  to  the  tilt  axis)  of  an  atomically 
relaxed  structure  of  the  2  =  5{310}<001)  symmetric  grain 
boundary  of  Si  in  its  ground  state  (GBl).  In  the  actual  3D 
structure,  the  atoms  in  the  bulk  lie  on  four  different  planes.  The 
letters  a-g  denote  various  sites  at  which  segregation  of  isolated 
As  atoms  and  As  dimers  are  investigated.  The  dislocation  cores 
comprising  this  grain  boundary  are  of  the  pure  edge  type. 

boundaries.  Figures  1  and  2  show  labels  for  the  sites 
where  As  atoms  were  placed:  [a]  through  [e]  are  sites  in 
the  GBl  grain  boundary,  whereas  [f]  and  [g]  are  sites  in 
the  bulk;  [a']  through  [f']  denote  corresponding  sites  for 
GB2.  Symmetry  has  been  used  to  reduce  the  number  of 
possible  distinct  As  sites  and  site  pairs.  Thus  using  the 
symbol  to  indicate  symmetry  equivalence,  we  have 
for  GBl:  a  ~  c  and  d  e.  It  follows  that  for  GBl 
the  possible  distinct  sites  for  atom  segregation  are  [a] 
([c]),  [b],  and  [d]  ([e]).  The  distinct  site  pairs  for  dimer 
segregation  are  [a,c],  [a,  d],  [b,c],  and  [d,e].  In  addition, 
we  have  also  studied  the  dimer  [f,g]  where  two  atoms 
are  placed  at  nearest-neighbor  sites  in  the  bulk  crystal. 
For  GB2  we  have  the  equivalence  fc'  ~  c'  e'.  The 
resulting  distinct  site  pairs  are  [a',d'],  [d',e'],  [b^c'],  and 
[a',c']. 

The  results  for  isolated  As  substitutionals  are  shown  in 
the  top  halves  of  Tables  I  and  IL  We  find  that  all  sites 
except  [a]  ([c])  on  GBl  and  [a^]  on  GB2  have  a  binding 
energy  of  '^O.l  eV,  the  same  value  obtained  for  isolated 


FIG.  3.  Grain  boundary  of  Fig.  2  in  a  metastable  state  (GB2) 
in  the  same  projection  view.  The  letters  a'-f'  denote  various 
sites  at  which  segregation  of  isolated  As  atoms  and  As  dimers 
are  investigated.  The  dislocation  cores  comprising  this  grain 
boundary  have  mixed  screw  and  edge  characters. 
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TABLE  I.  Segregation  energies  (normalized  per  As  atom)  for  isolated  As  atoms  and  As 
dimers  placed  at  various  sites  on  GB 1 .  The  equilibrium  As- As  distance  for  the  dimers  is  also 
indicated.  The  sites  are  described  in  Fig.  2. 


As  site(s) 

Segregation  energy 
eV/As  atom 

As-As  distance  (A) 
(dimer  geometries) 

As-As*  distance  (A) 
(dimer  geometries) 

[f] 

0.00 

[a],[c] 

0.03 

[b] 

0.11 

[d],  [e] 

0.11 

[f,g] 

-0.01 

2.71 

4.59 

[a,d] 

0.10 

2.71 

5.20 

[d,el 

0.19 

2.79 

4.70 

[a,c] 

0.20 

2.89 

3.76 

[b,c] 

0.32 

3.43 

3.43 

As  atoms  in  a  similar  boundary  in  Ge  [8].  The  smaller 
binding  of  [a]  can  be  explained  from  the  similarity  of 
its  environment,  characterized  by  the  surrounding  bond- 
length  and  bond-angle  distribution,  to  that  of  a  bulk  site 
[f]  or  [g],  while  an  opposite  effect  occurs  for  [a'].  In  all 
cases  of  isolated  As  substitutionals,  the  lattice  is  found  to 
undergo  only  a  very  small  relaxation,  in  agreement  with 
the  results  of  Ref.  [8]. 

We  now  turn  to  the  dimer  configurations.  As  we 
noted  earlier,  we  performed  calculations  with  both  mini¬ 
mal  and  double  periodicity  in  the  grain-boundary  plane, 
corresponding  to  dimer  chains  and  isolated  dimers,  re¬ 
spectively.  The  latter  calculations  are  extremely  time 
consuming  even  on  the  Paragon  supercomputer  and 
were  therefore  performed  only  for  selected  pairs  of  sites. 
We  found  that  dimer  formation  in  the  grain  boundary  is 
energetically  favored.  If  two  As  atoms  are  placed  at  neigh¬ 
boring  substitutional  sites  in  the  bulk  crystal,  they  repel 
each  other  seeking  to  achieve  threefold  coordination 
(pair  [f,e]  in  Table  I).  The  equilibrium  As-As  distance 
is  2.71  A  compared  to  a  Si-Si  bond  length  of  2.35  A. 
The  overall  energy  goes  up  by  a  tiny  amount  (0.01  eV 
per  As  atom)  as  compared  with  isolated  substitutional 
atoms  because  of  the  elastic  energy  cost.  In  contrast. 
As  dimers  in  the  grain  boundary  lead  to  an  overall 
lowering  of  the  energy.  The  selected  calculations  we 
performed  for  the  isolated  dimers  yielded  net  binding  of 
order  0.05-0.2  eV/atom,  indicating  that  the  elastic  energy 
cost  in  the  grain  boundary  can  be  smaller  than  in  the 


bulk.  Dimer  formation  in  the  grain  boundary  is  the  result 
of  repulsion  between  neighboring  As  atoms  and  occurs 
because  this  repulsion  can  be  accommodated  easier  in  the 
grain  boundary  than  in  the  bulk  crystal. 

The  results  for  chains  of  As  dimers  are  even  more 
dramatic  and  are  displayed  in  detail  in  Tables  I  and  n  for 
the  two  grain-boundary  structures,  respectively.  We  see 
that  segregation  energies  range  from  0.1  to  0.5  eV/atom, 
the  latter  being  in  agreement  with  the  measured  values 
[1-3].  In  Tables  I  and  II,  the  third  column  contains  the 
As-As  distance  in  the  dimer,  which  is  to  be  compared 
with  the  normal  Si-Si  distance  of  2.35  A.  The  fourth 
column  contains  the  As-As*  distance  between  As  atoms 
of  neighboring  dimers  (dimers  in  neighboring  supercells). 

We  note  three  classes  of  results:  (i)  Cases  where  the 
As-As*  distance  is  significantly  larger  than  the  As-As 
distance,  suggesting  that  the  dimers  in  the  chain  are  fairly 
well  separated.  The  segregation  energy  is  small,  less  than 
0.2  eV/atom,  comparable  to  that  of  truly  isolated  dimers 
that  we  discussed  earlier,  (ii)  Cases  where  the  As-As  and 
As-As*  distances  are  comparable  but  different  (e.g.,  [a,  c] 
in  GBl  and  [a',c']  in  GB2)  where  the  segregation  energy 
ranges  from  0.2  to  0.5  eV/atom.  (iii)  A  case  where  the 
As-As  and  As-As*  distances  are  identical  ([b,c]  in  GBl), 
corresponding  to  a  fully  ordered  chain  of  As  atoms, 
with  an  intermediate  segregation  energy  of  0.32  eV/atom. 
Figure  4  displays  the  electronic  charge  density  in  a  slice 
passing  through  the  plane  containing  the  As  atoms  in  the 
[b, c]  geometry  of  GBl.  There  is  no  significant  charge 


TABLE  II.  Segregation  energies  and  As-As  distance  at  various  sites  on  GB2  (Fig.  3). 


As  site(s) 

Segregation  energy 
eV/As  atom 

As-As  distance  (A) 
(dimer  geometries) 

As-As*  distance  (A) 
(dimer  geometries) 

[f'] 

0.00 

[ba[c'].[e'] 

0.12 

[d'] 

0.13 

[a'] 

0.22 

[a'.d'] 

0.08 

2.43 

5.25 

[d'.e'] 

0.09 

2.72 

4.84 

[b',c'] 

0.11 

2.42 

4.29 

[a',c'] 

0.52 

2.76 

3.54 
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FIG.  4(color).  A  2D  plot  of  the  electronic  charge  density  on  a 
slice  passing  through  the  perfect  As  chain  [b,c]  in  GBl.  The 
As  atoms  are  magenta  and  Si  atoms  gray.  The  color  scheme  is 
red  yellow  — ^  green  — ►  light  blue — ►  deep  blue  for  maximum 
to  minimum  charge  density. 


in  the  region  between  two  neighboring  As  atoms  in  the 
chain,  proving  conclusively  the  threefold  coordination  of 
each  As  atom.  The  charge  density  for  all  other  dimer 
configurations  is  qualitatively  similar. 

An  examination  of  the  local  three-dimensional  topolo¬ 
gies  corresponding  to  Fig.  2  and  3  suggests  that  the  most 
stable  dimer  geometries  in  a  chain  are  the  ones  in  which 
one  or  both  of  the  component  As  atoms  can  relax  into 
the  dislocation  cores,  where  the  average  atomic  density  is 
lower  than  the  crystalline  bulk.  Such  relaxation  reduces 
the  strain  in  the  Si  backbonds,  thereby  yielding  a  larger 
segregation  energy.  It  is  interesting  to  note  that  the  most 
stable  chain  of  dimers  [a',  c']  in  GB2  has  its  component 
atoms  on  two  different  dislocation  cores,  a  configuration 
possible  only  in  relatively  large-angle  grain  boundaries, 
while  the  chain  of  dimers  [b,c]  on  GBl  has  only  one  of 
its  atoms  (c)  lying  on  a  dislocation  core,  and  could  occur 
in  any  dislocation  core  of  the  perfect  edge  type. 

Finally,  it  is  particularly  interesting  to  note  that  the 
highest  segregation  energy  in  GB2  is  significantly  larger 
than  in  GBl  and  might  actually  drive  a  structural  transfor¬ 
mation.  In  other  words,  the  injection  of  a  high  concen¬ 
tration  of  As  into  polycrystalline  Si  may  convert  a  grain 
boundary  with  perfect  edge  dislocations  into  one  with 
mixed  dislocations.  Similar  solute-induced  grain  bound¬ 
ary  transformations  have  long  been  known  to  occur  in 
metals  [18],  but  we  are  not  aware  of  any  prior  reports 


in  semiconductors.  From  the  relative  total  energies  of  the 
two  grain  boundaries  we  estimate  that  the  transformation 
from  GBl  to  GB2  would  require  a  critical  As  concentra¬ 
tion  of  —19%  in  the  column  of  the  favored  dimer  sites. 
However,  such  a  transformation  involves  a  relative  shift 
(sliding)  of  the  two  grains  at  the  boundary,  and  the  above 
estimate  does  not  take  into  account  any  elastic  energy  cost 
that  may  be  required  to  maintain  integrity  at  triple  junc¬ 
tions  during  the  sliding  process. 

In  summary,  cooperative  phenomena  involving  chains 
of  threefold-coordinated  As  atoms  or  dimers  result  in  much 
larger  segregation  energies  than  isolated  As  substitutionals. 
Segregation  energies  thus  obtained  are  in  agreement  with 
experimental  values.  This  provides  a  mechanism  for  As 
segregation  that  does  not  require  the  presence  of  steps  or 
other  defects.  Chains  of  As  dimers  in  mixed  dislocation 
cores  have  lower  energies,  raising  the  intriguing  possibility 
that  As  segregation  may  drive  a  structural  transformation 
in  grain  boundaries  containing  pure  edge  dislocations. 
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We  demonstrate  by  ab  initio  calculations  that  segregation  of  As  in  a  dislocation  core  in  Si  occurs  in 
the  form  of  an  ordered  chain  of  As  atoms  running  along  the  dislocation  pipe.  All  As  atoms  in  the 
chain  achieve  threefold  coordination  and  the  segregation  energy  is  close  to  1  eV  per  As  atom. 
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Plastically  deformed  materials  release  excess  strain  by 
creating  extended  defects  such  as  dislocations.  In  semicon¬ 
ductor  devices,  dislocations  can  severely  affect  the  behavior 
of  dopant  impurities  because:  (1)  the  dislocation  cores  may 
provide  a  fast  diffusion  pathway  to  impurities,  which  se¬ 
verely  influences  the  dopant  profile  in  an  uncontrolled  way, 
and  (2)  impurities  may  get  trapped  in  the  core  regions  lead¬ 
ing  to  preferential  segregation  and  electrical  deactivation  of 
the  impurity. 

It  has  generally  been  believed  that  impurities  get  trapped 
at  sites  with  imperfect  bond  order  in  dislocation  cores.  How¬ 
ever,  experimental and  theoretical"^"'^  work  has  estab¬ 
lished  that  dislocations  and  grain  boundaries  in  Si  and  Ge 
actually  reconstruct  in  such  a  manner  that  all  atoms  are  four¬ 
fold  coordinated.  Therefore,  a  single  impurity  atom  occupy¬ 
ing  a  substitutional  site  in  the  core  does  not  experience  any 
appreciable  difference  from  its  environment  in  the  bulk  ma¬ 
terial  away  from  the  core,  resulting  in  small  segregation  en¬ 
ergies.  For  example,  in  the  case  of  As  segregation  in  a  sym¬ 
metric  tilt  boundary  in  Ge,  the  segregation  energy  was  found 
to  be  only  of  order  0.1  eV.'^  However,  in  a  very  recent  ab 
initio  study we  found  that  As  impurities  form  ordered 
chains  of  As  dimers  (a  periodically  repeated  array  of  two 
substitutional  As  atoms  at  neighboring  sites),  or  in  some 
cases,  perfectly  ordered  chains  of  As  atoms  along  the  grain 
boundary  cores,  with  segregation  energies  as  high  as  0.5  eV 
per  As  atom.  In  such  chain  configurations,  each  As  atom 
relaxes  away  from  its  partner  atom  in  the  dimer  through  re¬ 
pulsion,  and  attains  its  preferred  threefold  coordination, 
thereby  lowering  energy.  It  was  also  found  that  isolated  As 
dimers  lead  to  threefold  coordination,  but  the  excess  binding 
energy  over  single  substitutional  As  atoms  is  small,  ~0.1  eV 
per  atom,  and  sometimes  there  is  no  energy  gain  at  all.  It 
follows  that  the  energy  gain  from  attaining  threefold  coordi¬ 
nation  is  largely  cancelled  by  the  elastic  energy  cost  of  dis¬ 
torting  the  Si  backbonds.  However,  a  significant  additional 
energy  gain  occurs  through  chain  formation. 

In  this  letter,  we  show  that  As  segregation  at  an  isolated 
dislocation  core  in  Si  occurs  in  the  same  manner  as  in  grain 
boundaries,  namely  in  the  form  of  ordered  chains  of 
threefold-coordinated  As  atoms.  In  fact,  the  strains  associ¬ 
ated  with  isolated  dislocation  cores  are  in  general  much 
larger  than  in  grain  boundaries.  These  strains  lead  to  larger 
segregation  energies  in  a  dislocation  core  than  in  the  grain 
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boundary,  as  high  as  ^--0.9  eV  per  As  atom.  In  earlier  work 
using  cluster  calculations,  Jones  et  al?  found  that  As  and  P 
dimers  in  a  dislocation  core  gain  substantial  energy  by 
achieving  threefold  coordination.  However,  chains  of 
threefold-coordinated  impurities  were  not  considered. 

Before  delving  into  computational  details,  a  brief  intro¬ 
duction  to  dislocation  terminology  is  in  order.  The  major 
operative  slip  system  in  diamond  cubic  structures  involves 
the  glide  of  {111}  planes  along  (110)  directions.  Well  sepa¬ 
rated  perfect  glide  dislocations  lie  primarily  along  (110)  di¬ 
rections,  and  are  either  screw  or  60°  dislocations,  the  latter 
deriving  its  name  from  the  fact  that  the  angle  between  the 
dislocation  line  direction  and  the  Burgers  vector  b  is  60°.  In 
materials  with  low  stacking  fault  energy,  the  60°  dislocation 
is  often  found  to  be  dissociated  into  two  Shockley  partial 
dislocations  separated  by  a  stacking  fault.  These  two  Shock- 
ley  partials  have  their  Burgers  vectors  (b— Kll2))  aligned, 
respectively,  at  30°  and  90°  to  their  (110)  line  direction,  and 
are  commonly  referred  to  as  the  30°  and  90°  partials.'^  In 
this  work,  we  have  chosen  the  90°  partial  dislocation  core  in 
Si  as  a  concrete  system  to  perform  our  segregation  study. 

As  in  our  previous  work,'^  the  total  energy  calculations 
and  structural  relaxations  were  carried  out  using  density 
functional  theory  with  the  exchange  and  correlation  energy 
treated  in  the  local  density  approximation.  We  used  a  non- 
cubic  periodic  supercell  with  two  oppositely  oriented  90° 
partial  cores  separated  by  a  distance  of  13  A  (Fig.  1).  An 
energy  cutoff  of  150  eV  was  used,  and  the  Brillouin  zone 
integration  was  performed  using  two  special  k  points,  chosen 
according  to  the  Monkhorst-Pack  scheme.'^  Atoms  were  re- 


FIG.  I.  View  along  the  dislocation  line  of  a  periodic  supercell  containing 
two  oppositely  oriented  90°  partials  with  the  asymmetric  reconstmction;  the 
supercell  for  the  symmetric  reconstruction  appears  almost  identical  in  the 
same  view.  The  supercell  contains  64  atoms  and  the  cores  are  separated  by 
13  A. 
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FIG.  2.  Bonding  structure  of  the  90°  partial  dislocation  dipole  in  the  slip 
plane,  shown  by  doubling  the  supercell  periodicity  in  Fig.  1  along  the  dis¬ 
location  line  for  :  (a)  asymmetric  reconstruction  and  (b)  symmetric  recon¬ 
struction.  In  the  three-dimensional  geometry,  the  asymmetric  reconstruction 
[Fig.  2(a)]  has  all  atoms  fourfold  coordinated,  while  the  symmetric  recon¬ 
struction  [Fig.  2(b)]  has  some  atoms,  distinguished  by  lighter  shading,  quasi¬ 
fivefold-coordinated.  See  text. 

laxed  in  small  steps  until  the  magnitude  of  the  largest  force 
component  was  less  than  0,07  eV/A. 

Previous  studies  found  two  models  of  reconstruction  in 
the  90®  partial  dislocation  core:  (1)  the  asymmetric  recon¬ 
struction,  in  which  the  mirror  symmetry  along  the  dislocation 
line  is  broken/^’^^  and  (2)  the  symmetric  reconstruction,  in 
which  the  mirror  symmetry  is  kept  intact.  The  asymmetric 
structure  has  all  the  atoms  fourfold  coordinated,  while  the 
symmetric  structure  has  two  quasi-fivefold-coordinated  at¬ 
oms  per  periodic  segment  of  the  core.  Both  these  recon¬ 
structed  cores  look  almost  the  same  in  the  projection  along 
the  dislocation  line  (Fig.  1),  but  are  clearly  distinguishable  in 
the  bonding  structures  in  the  slip  plane  (Fig.  2).  From  ab 
initio  relaxations  we  found  that  the  asymmetric  structure  was 
stable,  while  the  symmetric  structure  spontaneously  trans¬ 
formed  into  the  asymmetric  structure,  in  agreement  with  an 
earlier  ab  initio  study.  All  our  segregation  studies  were 
thus  performed  exclusively  on  the  asymmetric  structure. 

Single  As  atoms  and  As  dimers  were,  respectively, 
placed  at  various  sites  and  nearest  neighbor  site  pairs  of  the 
dislocation  core,  as  indicated  in  Fig.  3.  In  all  cases,  the  en¬ 
ergy  of  the  relaxed  structure  was  calculated  relative  to  the 
geometry  in  which  the  As  atoms  in  the  core  were  exchanged 
with  Si  atoms  in  the  bulk,  i.e.,  away  from  the  core.  The 
negative  of  this  energy,  known  as  the  segregation  energy,  is 


FIG.  3.  Various  sites  of  an  isolated  90°  partial  core,  viewed  along  the 
dislocation  line,  on  which  single  As  atoms  and  As  dimers  are  placed.  Seg¬ 
regation  energies  are  listed  in  Table  I. 


TABLE  I.  Segregation  energies  of  single  As  atoms  and  As  dimers  placed 
substitutionally  at  various  sites  of  an  asymmetrically  reconstructed  core.  For 
the  dimer  geometries  the  relaxed  distance  between  the  two  dimer  companion 
atoms  and  the  corresponding  distance  between  the  first  As  atom 

and  the  periodic  image  of  its  dimer  companion  (^as-as*)  are  also  indicated. 
See  Fig.  3  for  site  denominations. 


As  site  (s) 

Segregation  energy 
(eV/As  atom) 

4as-as  (A) 
(Dimer  geometries) 

4 As- As*  (^) 

(Dimer  geometries) 

Bulk 

0.00 

[d] 

-0.15 

[a] 

0.14 

[b],  [c] 

0.33 

[d.f] 

0.01 

2.57 

2.42 

[d.e] 

0.13 

2.90 

4.53 

[b,c] 

0.88 

2.87 

2.95 

listed  in  Table  I  for  all  different  As  configurations  consid¬ 
ered. 

Let  us  first  consider  the  single  As  atoms  in  the  disloca¬ 
tion  core.  Four  different  sites  were  chosen.  The  strain  distri¬ 
bution  is  very  different  at  different  sites,  as  is  apparent  from 
comparing  the  site-associated  bond  lengths  with  the  bulk 
Si-Si  bond  length  of  2.35  A,  making  the  segregation  energy 
strongly  site  dependent.  Thus  the  strain  is:  (i)  compressive  at 
site  a,  with  bond  lengths  (2.29,  2.31,  2.31,  and  2.34)  A;  (ii) 
tensile  at  site  d,  with  bond  lengths  (2.37,  2.41,  2.43,  and 
2.44)  A;  (iii)  mixed  compressive/tensile  at  sites  b  and  c  with 
bond  lengths  (2.31,  2.34,  2.41,  and  2.43)  A.  The  resulting 
segregation  energy  for  As  is  largest  for  sites  b  and  c,  mod¬ 
erate  but  positive  at  site  a,  and  even  negative  for  site  d.  The 
largest  segregation  energy  of  0.33  eV  (sites  b,c)  is  much 
larger  than  the  average  binding  at  a  grain  boundary. 
Also,  the  wide  variation  of  segregation  energy  from  site  to 
site  is  to  be  contrasted  with  a  much  smaller  dispersion  in  a 
grain  boundary,  where  the  strain  distribution  is  much  more 
uniform.  In  all  cases  of  single  substitutional  As,  the  relax¬ 
ation  from  the  initial  Si  site  is  very  small,  just  as  in  a  grain 
boundary. 

Segregation  studies  of  periodic  chains  of  As  dimers  yield 
more  dramatic  results.  To  interpret  these  results,  it  should  be 
reemphasized  that  two  As  atoms  placed  on  nearest  neighbor 
Si  sites,  which  we  call  As  dimers,  seek  to  repel  and  relax 
away  from  each  other.  In  this  way,  each  As  atom  achieves  its 
preferred  threefold  coordination,  leaving  all  the  Si  atoms 
fourfold  coordinated.^^  The  three  largest  bonds  in  the  dislo¬ 
cation  core,  i.e.,  [d,f],  [d,e],  and  [b,c]  were  chosen  for  inves¬ 
tigation.  Table  I  displays  the  segregation  energies  for  these 
three  dimers.  It  also  lists,  for  each  relaxed  geometry,  the 
distance  of  an  As  atom  from:  (a)  its  dimer  companion  (col¬ 
umn  3),  and  (b)  the  nearest  periodic  image  (As*)  of  the 
dimer  companion  in  the  dislocation  line  direction  (column 
4).  The  distinct  behaviors  of  the  three  dimers  are  clearly 
evident  from  Table  I,  as  discussed  below: 

(i)  [d,f]:  In  this  case  the  orientation  of  the  dimer  is  such 

that  site  d  is  the  nearest  neighbor  of  both  f  and  its 
nearest  periodic  image  /*  in  the  dislocation  line  di¬ 
rection.  Thus,  in  the  periodic  chain  of  dimers,  all  the 
As  atoms  are  too  close  to  each  other.  This  allows  only 
a  small  stretching  of  the  As-As  separation,  to  only 
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2.57  A,  while  the  As-As*  distance  (2.42  A)  changes 
very  little  from  the  equilibrium  Si-Si  bond  distance 
of  2.35  A.  Consequently,  the  relaxed  structure  con¬ 
sists  of  a  chain  of  essentially  fourfold-coordinated  As 
atoms,  resulting  in  negligible  binding. 

(ii)  [d,e]:  In  this  case  the  orientation  of  the  dimer  is  nearly 
perpendicular  to  the  dislocation  line  direction,  and 
sites  d  and  e*  are  only  second  neighbors.  This  allows 
the  As  dimer  atoms  to  relax  away  from  each  other, 
stretching  the  As-As  separation  to  2.90  A,  thereby 
leading  to  threefold  coordination  of  each  As  atom. 
However,  the  transverse  orientation  of  the  dimer 
keeps  it  well  separated  from  its  periodic  images  in  the 
dislocation  line  direction,  as  is  evident  from  the  large 
As-As*  distance.  We  therefore  have  an  array  of 
nearly  isolated  As  dimers.  This  leads  to  a  positive  but 
low  segregation  energy  (0.13  eV  per  As)  ,  just  as  in  a 
grain  boundary. 

(iii)  [b,c]:  In  this  case  the  dimer  is  oriented  such  that  al¬ 
though  b  and  c*  are  second  nearest  neighbors,  their 
distance  before  relaxation  is  much  less  than  the  sepa¬ 
ration  of  d  and  e*  in  case  (ii).  Consequently,  the 
As-As  stretch  is  accompanied  by  As  and  As*  getting 
closer  to  each  other.  In  the  final  relaxed  geometry,  the 
As-As*  separation  (2.95  A)  is  only  slightly  larger 
than  the  As-As  separation  (2.87  A)  and  we  have  a 
nearly  perfect  chain  of  As  atoms  leading  to  a  large 
segregation  energy  (0.88  eV  per  As). 

In  order  to  elucidate  further  the  source  of  energy  gain 
through  the  formation  of  dimer  chains,  we  have  examined 
the  electronic  energy  levels  in  the  As  dimers  in  the  disloca¬ 
tion  core  and  compared  them  with  those  of  isolated  As  atoms 
in  the  bulk.  As  one  would  expect,  the  isolated  As  atom  has  a 
shallow  donor  level  at  less  than  0. 1  eV  below  the  conduction 

band  edge.  In  the  [b,c]  dimer  chains,  on  the  other  hand,  each 
^  ■  20 
As  atom  has  an  electron  at  a  level  in  the  midgap  region. 

Further  analysis  indicates  that  the  level  in  the  gap  is  a  mem¬ 
ber  of  the  lone-pair  states  that  As  atoms  have  when  they  are 
threefold  coordinated.  Because  of  the  close  proximity  of  the 
threefold  coordinated  As  atoms,  these  lone-pair  states  split, 
with  half  of  them  in  the  band  gap  and  half  of  them  in  the 
valence  band.  Thus,  the  gain  in  energy  could  be  viewed  as  a 
result  of  the  shallow  donor  level  being  driven  deep  into  the 
band  gap  by  the  lattice  relaxation  accompanying  dimer  for¬ 
mation,  or,  equivalently,  as  arising  from  the  fact  that  dimer 
formation  leads  to  lone-pair  states  that  are  lower  in  energy 
than  the  states  that  are  available  when  an  As  atom  is  fourfold 
coordinated. 


In  summary,  using  the  90°  partial  dislocation  as  a  con¬ 
crete  example,  we  find  that  As  likes  to  decorate  a  core  in  Si 
in  an  ordered  chain  fashion.  All  As  atoms  in  the  chain 
achieve  threefold  coordination  through  repulsion  between  al¬ 
ternating  As  atoms.  Overall  behavior  of  single  impurities, 
isolated  dimers,  and  chains  in  the  dislocation  core  is  similar 
to  the  behavior  we  found  in  a  Si  grain  boundary.  However, 
the  segregation  energy  of  both  single  As  atoms  and  those  in 
a  chain  are  much  larger  in  the  dislocation  core,  with  the 
result  for  the  chain  being  as  high  as  0.9  eV  per  As  atom. 
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Appendix  D 

Coupled-Barrier  Diffusion:  The  Case  of  Oxygen  in  Silicon 

Madhavan  Ramamoorthy  and  Sokrates  T.  Pantelides 
Department  of  Physics  and  Astronomy,  Vanderbilt  University,  Nashville,  Tennessee  37209 

(Received  9  August  1995) 

Oxygen  migration  in  silicon  corresponds  to  an  apparently  simple  jump  between  neighboring  bridge 
sites.  Yet  extensive  theoretical  calculations  have  so  far  produced  conflicting  results  and  have  failed  to 
provide  a  satisfactory  account  of  the  observed  2.5  eV  activation  energy.  We  report  a  comprehensive 
set  of  first-principles  calculations  that  demonstrate  that  the  seemingly  simple  oxygen  jump  is  actually  a 
complex  process  involving  coupled  barriers  and  can  be  properly  described  quantitatively  in  terms  of  an 
energy  hypersurface  with  a  “saddle  ridge”  and  an  activation  energy  of  ~2.5  eV.  Earlier  calculations 
correspond  to  different  points  or  lines  on  this  hypersurface. 


PACS  numbers:  66.30.Jt,  31.15.Ar,  61.72.-y,  81.60.Cp 

Oxygen  in  silicon  has  long  been  known  to  occupy 
a  bridge  position  between  neighboring  Si  atoms,  with 
an  Si-O-Si  configuration  similar  to  those  in  Si02  [1,2]. 
Its  diffusion,  measured  to  have  an  activation  energy  of 
2.5  eV  [3],  is  generally  believed  to  consist  of  simple 
jumps  between  neighboring  bridge  positions  on  the  (110) 
plane  defined  by  the  corresponding  Si-Si  bonds  (Fig.  1). 
In  terms  of  the  angle  do  defined  in  Fig.  1 ,  the  midpoint  of 
the  jump  is  at  Oq  =  90°. 

Most  calculations  to  date  [4—8]  assumed  such  a  simple 
adiabatic  jump,  with  reflection  symmetry  about  the  ver¬ 
tical  axis  shown  in  Fig.  1.  Thus,  the  saddle  point  was 
assumed  to  have  the  O  atom  at  Oq  =  90°  and  the  central 
Si  atom  at  ^si  =  90°.  The  remaining  degrees  of  free¬ 
dom  and  the  positions  of  the  other  Si  atoms  were  deter¬ 
mined  by  total-energy  minimization.  The  resulting  total 
energy,  measured  from  the  energy  of  the  equilibrium  con¬ 
figuration,  represents  the  adiabatic  activation  energy  for 
diffusion.  Some  authors  [4,5]  reported  activation  ener¬ 
gies  around  2.5  eV,  while  others  [6—8]  reported  smaller 
values  ranging  from  1.2  to  2.0  eV. 

In  Ref.  [8],  Needels  et  al.  found  a  value  of  1.8  eV  and 
attributed  the  discrepancy  with  experiment  to  dynamical 
phenomena,  i.e.,  the  neighboring  Si  atoms  do  not  relax 
fully  along  the  O  trajectory.  TTiey  reported  model  dy¬ 
namical  calculations  for  a  “generic”  nonadiabatic  path  in 
which  the  O  atom  was  given  an  initial  “kick,”  i.e.,  an  ini¬ 
tial  velocity  corresponding  to  a  kinetic  energy  of  2.0,  2.3, 
or  2.7  eV.  They  found  that  when  the  kick  energy  was 
<2.5  eV,  the  O  atom  went  past  the  saddle  point  but  then 
returned  to  the  original  bridge  position.  When  the  kick 
was  >2.5  eV  the  O  atom  migrated  to  the  next  bridge  site. 
They  concluded  that  their  results  suggested  that  dynamical 
effects  are  important  in  O  migration,  but  did  not  constitute 
definitive  evidence. 

In  a  recent  paper,  Jiang  and  Brown  (JB)  [9]  sought  to 
resolve  the  issue  by  exploring  the  entire  migration  path. 
They  performed  total-energy  minimizations  by  stepping 
the  oxygen  atom  from  one  bridge  site  to  the  next.  They 
found  that  the  total  energy  attains  a  value  of  only  —1.2  eV 
at  —  90°,  but  then  keeps  rising  to  a  maximum  value 

003 1-9007/ 96/76(2)/267(4)$06.00 


(saddle  point)  of  -2.5  eV  at  =  113°.  In  addition, 
they  computed  the  diffusion  constant  and  found  it  to 
agree  very  well  with  experiment  over  12  decades.  They 
concluded  that  the  saddle  point  of  O  migration  is  past  the 
midpoint  of  the  path  and  that  their  results  account  for  all 
the  experimental  data. 

At  first  glance,  JB’s  results  account  nicely  for  the  ex¬ 
perimental  data  without  the  need  to  invoke  dynamical  ef¬ 
fects.  Nevertheless,  the  pronounced  asymmetry  of  JB’s 
total-energy  profile  about  do  =  90°  raises  a  serious  ques¬ 
tion:  If  the  global  minimum  of  the  total  energy  was 
indeed  obtained  at  each  point  of  the  O  path,  the  energy 
profile  would  be  symmetric  about  do  =  90°.  Clearly,  JB’s 
minimization  procedure  yielded  a  secondary  minimum  for 
each  do  >  90°,  not  the  global  minimum.  The  energy  at 
the  global  minimum  for  90*  +  a  is  by  symmetry  equal  to 
that  at  90°  —  a.  If  an  energy  profile  were  constructed  us¬ 
ing  global  minima  along  the  entire  path,  it  would  have  a 
maximum  of  only  1.2  eV  at  do  =  90°.  This  value  would 
be  in  poor  agreement  with  experiment.  We  conclude  that 
there  is  still  no  satisfactory  account  of  the  observed  2.5  eV 


I 


I 

I 

I 

I 

FIG.  1.  The  geometry  of  O  migration  in  a  (110)  plane.  The 
solid  dots  are  the  nominal  positions  of  Si  atoms  in  the  perfect 
crystal.  The  open  circles  show  the  positions  of  the  Si  and  O 
atoms  in  the  equilibrium  configuration. 
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activation  energy,  or  of  the  mutually  conflicting  theoretical 
results  published  so  far  on  oxygen  diffusion  in  silicon. 

In  this  paper,  we  report  a  series  of  first-principles  total- 
energy  calculations  which  show  that  the  process  of  O  mi¬ 
gration  is  far  more  complex  than  has  been  recognized  so 
far,  but  is  still  adiabatic.  During  the  migration  process, 
both  the  O  atom  and  the  central  Si  atom  perform  jumps 
and  face  large  barriers.  As  a  result,  a  quantitative  de¬ 
scription  of  the  process  requires  a  calculation  of  the  total- 
energy  hypersurface  as  a  function  of  the  positions  of  both 
atoms.  We  will  show  slices  of  this  hypersurface  that  re¬ 
veal  a  “saddle  ridge”  in  multidimensional  space.  The 
migration  process  is  adiabatic  and  occurs  along  a  multiplic¬ 
ity  of  paths  over  this  ridge  with  a  predominant  barrier  of 
~2.5  eV.  Finally,  we  find  that  the  results  of  earlier  au¬ 
thors  correspond  to  different  points  or  lines  on  the  hyper- 
surface. 

We  performed  calculations  using  density  functional  the¬ 
ory  and  the  local-density  approximation  for  exchange  and 
correlation,  using  the  form  for  the  exchange-correlation  po¬ 
tential  given  by  Ceperley  and  Alder  [10].  The  ultrasoft 
pseudopotentials  of  Vanderbilt  [11]  were  used  for  Si  and 
O.  These  pseudopotentials  have  been  thoroughly  tested 
in  several  extensive  investigations  [12-14].  The  calcula¬ 
tions  employed  a  plane  wave  basis  set  and  converged  re¬ 
sults  were  obtained  with  an  energy  cutoff  of  25  Ry.  A  bcc 
supercell  with  32  Si  atoms  and  one  O  atom  was  used.  Each 
structure  was  relaxed  until  the  force  on  each  atom  was  less 
than  0.5  eV/A.  All  calculations  were  first  done  with  one 
special  k  point  at  (0.5, 0.5, 0.5)  in  the  irreducible  Brillouin 
zone  [15].  The  key  calculations  were  repeated  with  two  k 
points  at  (0.75, 0.25, 0.25)  and  (0.25, 0.25, 0.25)  [15].  The 
energy  differences  changed  at  most  by  about  0.2  eV,  with 
all  the  qualitative  results  obtained  with  one  k  point  being 
unchanged.  Hence,  the  results  with  one  k  point  were  taken 
to  be  converged  with  respect  to  i-point  sampling,  and  used 
in  all  the  figures  in  this  paper. 

Our  results  for  the  equilibrium  configuration  of  O, 
shown  schematically  in  Fig.  1,  are  in  agreement  with 
earlier  work  [2].  We  find  a  very  flat  minimum  at  Oq  ~ 
55  .  The  Si-O  bond  length  is  1.6  A,  the  Si-Si  length  is 
3.2  A  (compare  with  the  value  of  2.35  A  in  bulk  silicon), 
and  the  Si-O-Si  bond  angle  is  ~150°.  The  angle  6si  is 
also  —150°.  For  our  purposes  here,  the  key  point  is  that 
the  central  Si  atom  is  well  to  the  right  of  the  vertical 
symmetry  axis  (see  Fig.  1).  As  the  O  migrates  from  the 
left  bridge  position  to  the  one  on  the  right,  the  central  Si 
needs  to  move  from  the  right  to  the  left,  specifically  from 
^Si  ~  150°  to  ^si  ~  30°.  We  will  see  below  that  this 
swing  of  the  central  Si  atom  controls  the  dynamics  of  the 
oxygen  migration  because  the  Si  atom  has  to  overcome  a 
barrier. 

We  demonstrate  this  key  result  in  Fig.  2,  where  we 
plot  the  total  energy  of  the  system  as  a  function  of  ^si 
when  the  O  atom  at  Bq  =  90°  [16].  For  each  Osu  the 
total  energy  was  minimized  with  respect  to  Rq,  Rsi  (see 
Fig.  1),  and  the  positions  of  the  other  Si  atoms.  We 
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FIG.  2.  The  total-energy  variation  of  the  system,  as  a  function 
of  Os,  when  the  O  atom  is  at  0o  =  90°.  The  zero  of  the  energy 
in  this  and  in  subsequent  plots  is  taken  as  that  of  the  equilibrium 
configuration. 

see  that  it  costs  only  0.6  eV  to  place  the  O  atom  at  the 
midpoint  if  the  central  Si  atom  is  allowed  to  relax  freely 
either  to  the  left  or  to  the  right  [17].  As  we  saw  earlier, 
as  the  O  atom  moves  from  the  left  bridge  position  to 
the  one  on  the  right,  the  central  Si  atom  needs  to  swing 
from  the  tight  side  to  the  left  side.  Figure  2  shows  that, 
with  do  =  90°,  the  total  barrier  is  2.2  eV.  This  barrier 
corresponds  to  the  two  atoms  crossing  the  midpoints  of 
their  respective  paths  at  the  same  time.  It  could  be  argued 
that  this  configuration  constitutes  the  saddle  point,  as  was 
assumed  in  several  previous  investigations  [4-8].  The 
total  barrier  of  -2.2  eV  is  indeed  in  good  agreement 
with  the  experimental  value.  This  simple  result,  however, 
belies  an  enormous  complexity  which  we  unravel  below. 
The  basis  of  this  complexity  is  that  the  O  atom  and  the 
central  Si  atom  need  not  pass  through  the  midpoints  of 
their  paths  at  the  same  time. 

The  above  analysis  makes  it  clear  that  O  migration  needs 
to  be  described  in  at  least  a  two-dimensional  space  defined 
by  00  and  ^si  because  the  central  Si  atom  also  must  climb 
a  barrier.  This  barrier,  however,  is  not  simple,  but,  as 
shown  in  Fig.  2,  has  a  cusp  at  ^si  =  90°,  indicative  of 
a  Jahn-Teller-like  instability  with  two  symmetric  total- 
energy  manifolds.  These  two  manifolds  correspond  to 
the  central  Si  atom  being  to  the  left  or  the  right  of  the 
symmetry  axis,  being  bonded  to  the  respective  Si  atom  on 
the  left  or  the  right.  For  values  of  0o  other  than  90°,  the 
two  manifolds  are  not  symmetric.  In  Fig.  3  we  trace  the 
evolution  of  the  two  total-energy  manifolds  for  a  sequence 
of  values  starting  with  the  O  atom  near  its  equilibrium 
bridge  position  on  the  left  of  the  vertical  symmetry  axis 
(bottom  panel).  The  central  Si  atom  is  on  the  right  side 
of  the  axis  (the  lower-energy  manifold).  As  the  O  atom 
progresses  along  its  path  (higher  panels  in  Fig.  3),  the 
central  Si  atom  stays  in  the  right  manifold.  At  0o  =  90°, 
the  two  manifolds  cross  and  the  central  Si  atom  can  switch 
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FIG.  3.  A  series  of  plots  indicating  the  variation  of  the  total 
energy  as  a  function  of  ^si  for  several  values  of  ^o* 


manifolds  and  swing  over  to  the  left  of  the  axis,  so  that  both 
the  O  atom  and  the  central  Si  atom  can  head  for  their  final 
destinations.  The  total  barrier  for  this  process  is  2.2  eV. 

There  are  additional  possibilities,  however.  The  O 
atom  may  overshoot  the  midpoint  of  its  jump  without  the 
central  Si  atom  swinging  over.  The  relevant  total-energy 
manifolds  are  shown  in  the  upper  panels  of  Fig.  3.  The 
central  Si  atom  is  now  in  the  high-energy  manifold,  stuck 
on  the  “wrong”  side  of  the  vertical  axis.  Even  though 
the  manifolds  do  not  cross,  the  central  Si  atom  is  stable 
in  the  higher-energy  manifold  only  up  to  a  certain  value 
of  marked  by  the  solid  arrows.  At  those  points,  the 
calculations  show  that  the  central  Si  atom  collapses  to  the 
lower-energy  manifold,  i.e.,  swings  over  to  the  left  side 
of  the  axis.  No  matter  how  much  the  O  atom  overshoots 


(i.e.,  any  of  the  top  four  panels),  the  total  energy  needed 
for  the  central  Si  atom  to  swing  over  is  of  order  2.3- 
2.7  eV.  In  other  words,  the  O  atom  and  the  central  Si 
atom  need  not  move  in  concert  and  be  at  the  midpoints  of 
their  respective  paths  at  the  same  time.  They  can  move 
independently  and  still  face  a  net  barrier  of  ~2.5  eV. 

The  multiplicity  of  paths  is  best  illustrated  with  a  two- 
dimensional  plot  of  the  total  energy  as  a  function  of  and 
^Si .  For  clarity,  we  show  only  one  of  the  two  manifolds 
at  each  pair  (^o*  ^siX  namely,  the  one  corresponding  to 
coupled  adiabatic  migration.  The  surface  was  obtained  by 
interpolating  through  a  sizable  number  of  calculated  points. 
Note  the  flat  regions  corresponding  to  the  two  equilibrium 
configurations  and  the  steep  drops  that  correspond  to  the 
regions  of  the  solid  arrows  in  Fig.  3.  We  see  that  there 
is  a  “saddle  ridge”  with  a  net  energy  of  ^^2.5  eV  over  a 
considerable  range.  At  the  high  symmetry  point  on  the 
ridge  (^o  ==  90°,  ^si  =  90°),  the  total  energy  is  somewhat 
lower,  —2.2  eV,  but  this  smaller  value  corresponds  to  a 
small  fraction  of  all  possible  migration  paths  over  the  ridge 
(the  resolution  of  the  figure  is  limited  by  the  complexity  of 
the  surface  near  the  ridge).  There  are  two  classes  of  paths: 
those  in  which  the  O  atom  overshoots  the  midpoint  of  its 
path  with  the  central  Si  atom  trailing  and  those  in  which 
the  central  Si  atom  goes  over  the  midpoint  of  its  path  first 
with  the  O  atom  trailing.  Along  all  these  paths  the  distance 
between  the  O  atom  and  the  central  Si  atom  is  —1.7  A. 
Thus,  the  O-Si  bond  acts  like  a  pogo  stick  that  faces  a  net 
barrier  of  —2.5  eV  no  matter  how  it  turns  as  it  attempts  to 
change  its  tilt  from  the  left  to  the  right. 

The  collapse  from  one  manifold  to  the  other  indicated 
by  the  solid  arrows  in  Fig.  3  (steep  drops  in  Fig.  4)  was 
intriguing  enough  to  merit  further  investigation.  The  plots 
in  Figs.  3  and  4  were  constructed  by  picking  and  ^si 
and  then  letting  both  the  O  atom  and  the  central  Si  atom 
move  radially  until  the  energy  was  minimized.  It  turns 
out  that  the  two  manifolds  shown  in  Fig.  3  correspond  to 
two  fairly  distinct  regions  of  Rsi  values.  We  explored 
Rs\  values  between  these  two  regions  and  found  that  for 
each  6si  the  energy  as  a  function  of  Rs\  has  two  minima 
with  a  barrier  that  prohibits  the  central  Si  atom’s  motion 
from  one  minimum  to  another,  corresponding  to  a  switch 
between  the  two  energy  manifolds.  At  the  critical  ^si 
value  (solid  arrow  in  Fig.  3),  this  barrier  vanishes  and  the 
collapse  occurs.  The  evolution  of  this  radial  barrier  is  also 
quite  intriguing  and  will  be  discussed  further  in  a  longer 
article.  In  fact,  a  complete  description  of  O  migration 
requires  the  total  energy  as  a  function  of  four  coordinates 
on  the  (110)  plane:  (^o,  ^si»  ^Si)-  Figures  3  and  4 
represent  slices  through  this  hypersurface. 

We  now  turn  to  examine  the  earlier  theoretical  work  in 
light  of  the  present  work.  The  major  point  is  that  all  ear¬ 
lier  investigators  did  not  recognize  the  important  role  of 
the  central  Si  atom  in  the  migration  process.  Neverthe¬ 
less,  we  can  map  their  results  on  Fig.  4.  References  [4- 
8]  assumed  that  the  saddle  point  is  at  do  =  ^si  =  90°, 
shown  as  a  dot  in  Fig.  4.  Our  value  for  this  point  lies 
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FIG.  4.  The  total-energy  variation  as  a  function  of  do  and  ^si  • 
The  central  dot  and  the  hand-drawn  path  (in  the  direction  of  the 
arrowhead)  are  discussed  in  the  text. 

in  the  middle  of  the  range  of  earlier  theoretical  results 
(1.8 -2.5  eV)  [18].  This  spread  in  part  reflects  the  fact 
that  calculations  involving  oxygen  are  computationally 
extremely  demanding.  The  important  point  to  note  is  that 
in  Fig.  4,  the  paths  that  contain  the  point  do  =  ^si  =  90° 
constitute  a  small  portion  of  the  phase  space  of  all  paths 
going  over  the  ridge.  Most  of  the  ridge  is  somewhat 
higher,  ^2.5  eV,  close  to  the  observed  activation  energy. 

Figure  4  also  clarifies  the  “kick”  simulations  of  Needels 
€t  al.  [8]  and  JB  s  calculations.  Both  calculations  corre¬ 
spond  to  the  path  shown  by  the  hand-drawn  line  in  Fig.  4. 
Needels  et  al  gave  a  high  kinetic  energy  to  the  O  atom, 
whereas  JB  stepped  the  O  atom  gradually.  In  both  cases, 
the  O  atom  crossed  the  midpoint  of  its  path  before  the  cen¬ 
tral  Si  atom  did.  In  JB’s  case,  the  central  Si  atom  crossed 
oyer  when  ~  1 15°.  In  the  case  of  Needels  et  al,  for 
kick  energies  <2.5  eV,  the  O  atom  had  to  turn  back  be¬ 
cause  the  central  Si  atom  still  faced  a  barrier  and  could  not 
cross  over.  In  hindsight,  one  should  have  given  a  kick  to 
the  central  Si  atom.  In  any  case,  this  is  not  a  very  likely 
path  because  in  reality  both  the  O  atom  and  the  central 
Si  atom  are  vibrating  about  their  equilibrium  positions, 
attempting  to  overcome  their  respective  barriers.  As  the 
ridge  has  roughly  a  constant  height  over  a  considerable 
range,  any  point  at  an  energy  of  ~2.5  eV  may  be  fairly 
representative  of  the  entire  ridge.  This  might  explain  the 
good  agreement  obtained  by  JB  for  the  diffusion  constant 
with  experiment. 

Finally,  Needels  et  al  [8]  noticed  the  formation  of  a 
metastable  configuration  at  the  endpoint  of  the  O  path, 
when  ^e  kick  was  2.7  eV.  The  present  work  shows 
that  this  configuration  occurs  in  the  right-hand  manifold 
shown  in  the  upper  panels  of  Fig.  3.  The  metastable 
configuration  is  created  during  the  migration  process  for 
all  paths  that  do  not  cross  do  =  ^si  =  90°.  Contrary 
to  the  findings  of  Needels  et  al  [8],  this  configuration 
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collapses  to  the  equilibrium  configuration  at  the  end  of 
the  migration  process. 

In  summary,  we  have  shown  that  the  total-energy  varia¬ 
tion  dimng  O  migration  is  a  very  complex  hypersurface  in 
a  multidimensional  space.  We  have  shown  slices  of  this 
hypersurface  along  some  relevant  coordinates,  revealing 
seemingly  disconnected  manifolds.  We  believe  that  the 
calculations  we  have  done  so  far  have  captured  essentially 
all  of  die  very  complex  physics  of  the  seemingly  simple  O 
juinp  in  bulk  Si.  Our  results  suggest  that  such  complexity 
is  likely  to  be  present  whenever  migration  of  an  impurity 
involves  bond  breaking  and  rebonding  with  different  host 
atoms. 
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Abstract 


Experiments  have  established  two  modes  of  enhanced  oxygen  diffusion  in 
silicon,  one  of  unknown  origin  and  the  other  catalyzed  by  hydrogen.  Existing 
theories,  all  based  on  cluster  calculations,  have  yielded  contradictory  results 
for  both  modes.  We  report  state-of-the-art  first-principles  calculations  and 
determine  migration  pathways  for  a  self-enhanced  mode  via  oxygen  dimers 
and  for  the  H-enhanced  mode  via  oxygen-hydrogen  complexes.  The  concerted 
atomic  motions  are .  physically  transparent  and  the  corresponding  reduced 
activation  energies  agree  with  experiment. 
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Oxygen  diffusion  in  silicon  is  an  important  technological  process  because  it  leads  to  the 
formation  of  aggregates  (small  0  clusters,  Si02  precipitates)  during  the  thermal  processing 
of  crystalline  Si  wafers.^  The  published  experimental  data  exhibit  a  number  of  puzzles.^ 
While  the  diffusion  of  0  atoms  has  been  directly  observed^  to  take  place  with  an  activation 
energy  of  2.5  eV  over  a  wide  range  of  temperatures  (300-1200°C),  0  clustering  at  moderate 
temperatures  (350-450°  C)  takes  place  with  an  activation  energy  of  only  1.7  eV.”*”®  Moreover, 
the  observed  rates  of  clustering  at  these  temperatures  are  several  orders  of  magnitude 
higher  than  would  be  possible  by  means  of  the  diffusion  of  isolated  0  atoms.  To  resolve 
this  puzzle,  different  investigators  postulated  the  existence  of  different  candidates  for  a 
fast-diffusing  defect  complex:  the  O2  molecule,’'-®,  an  0-vacancy  complex,®  an  O-Si-self- 
interstitial  complex,^®  the  O  dimer.^^  No  definitive  experiment  or  theory  has  established  the 
complex’s  identity. 

Enhanced  0  diffusion  in  Si  has  also  been  shown  to  be  induced  by  H  in  a  separate  set  of 
investigations,^^  with  a  low  activation  energy  of  ~  1.2  eV.^®  The  catalytic  effect  of  H  has  been 
attributed  to  collisions  of  fast-diffusing  H  atoms  with  0  atoms^®  and  to  oxygen-hydrogen 
complexes.^'*  Again,  no  definitive  explanation  has  been  established. 

In  all  the  published  theoretical  work  on  enhanced  0  diffusion  in  Si,  the  infinite  Si  crystal 
was  approximated  by  a  small  cluster  containing  between  30  and  50  Si  atoms,  with  the 
impurity  atoms  placed  at  the  center.  Pathways  were  reported  for  “self-enhanced”  0  diffusion 
(i.e.  enhanced  0  migration  in  the  absence  of  any  other  impurity  or  native  defect)  via 
the  motion  of  0  dimers  and  for  H-enhanced  diffusion  via  the  motion  of  oxygen-hydrogen 
complexes.  The  results,  however,  are  mutually  contradictory.  For  self-enhanced  0  diffusion, 
Snyder  et  al.“,  using  the  semi-empirical  MINDO/3  method,  found  a  pathway  with  a  low 
activation  energy  of  ~  1.4  eV  via  the  migration  of  the  0  dimer.  However,  Ewels  et  al.^®,  using 
density  functional  theory,  did  not  find  a  lowered  activation  energy  along  the  same  pathway. 
For  H-enhanced  0  diffusion,  Estreicher^^  employed  the  approximate  PRDDO  method  and 
Jones  et  al.^®  employed  density  functional  theory  to  obtain  pathways  for  the  migration  of  the 
oxygen-hydrogen  complex  with  low  activation  energies  of  ~  1.5  eV.  However,  they  proposed 
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completely  different  configurations  for  the  equilibrium  structure  of  the  complex  and  different 
pathways  for  its  migration. 

In  this  paper  we  report  a  comprehensive  investigation  of  enhanced  migration  of  O 
in  Si  using  density  functional  theory,  the  local  density  approximation  for  the  exchange- 
correlation  potential  and  pseudopotentials.  The  defect  configurations  are  simulated  using 
supercells,  which  are  infinite  crystals  made  up  of  periodically  repeated  large  cells,  containing 
a  defect  at  the  center.  This  approach  has  been  established  over  the  last  decade  to  be 
reliable  for  addressing  problems  involving  the  energetics  and  diffusion  of  impurities  in  solids. 
In  particular,  a  few  years  ago,  this  method  was  adopted  by  one  of  the  present  authors 
and  his  then  coworkers  to  resolve  numerous  controversies  concerning  the  diffusion  of  H 
in  crystalline  Si.^^  More  recently,  the  present  authors  used  the  same  approach  to  resolve 
controversies  in  the  diffusion  of  the  0  atom  in  crystalline  Si.^®  Here,  we  report  the  results 
of  an  extensive  investigation  of  self-enhanced  and  H-enhanced  0  diffusion  in  Si.  We  show 
that  self-enhancement  occurs  with  a  low  activation  energy  of  about  1.5  eV,  through  the 
motion  of  O  dimers.  H-enhanced  diffusion  occurs  through  the  migration  of  oxygen-hydrogen 
complexes.  Our  results  for  the  equilibrium  configurations  for  the  positive,  neutral  and 
negative  charge  states  for  this  complex  bear  a  logical  similarity  to  the  published  results  for 
isolated  H  in  crystalline  Si.^’'  For  the  neutral  charge  state,  we  found  a  diffusion  pathway  with 
a  low  activation  energy  of  only  1.2  eV.  The  lowered  activation  energy  during  of  self-enhanced 
and  H-enhanced  0  diffusion,  in  coiriparison  with  the  migration  of  the  0  atom,  is  due  to  the 
saturation  of  a  dangling  bond  on  one  of  the  principal  Si  atoms  involved  in  the  migration 
process.  The  detailed  mechanisms  of  the  two  processes,  however,  are  very  different. 

The  theoretical  results  reported  in  this  article  were  obtained  using  density  functional 
theory  and  the  local-density  approximation  for  exchange  and  correlation.  The  ultra-soft 
pseudopotentials  of  Vanderbilt^^  were  used.  The  calculations  employed  a  plane  wave  basis 
set.  Converged  results  were  obtained  with  an  energy  cutoff  of  25  Ry.  A  32  site  bcc  supercell 
was  used,  with  one  k-point  at  (0.5, 0.5, 0.5)  in  the  irreducible  Brillouin  2one^°.  Each  structure 
was  relaxed  until  the  force  on  each  atom  was  less  than  0.5  eV /A. 
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In  the  following  paragraphs,  we  describe  migration  pathways  which  lead  to  enhanced  0 
migration  in  Si.  Given  the  complexity  of  the  energy  surface  describing  the  migration  of  each 
defect  complex,  we  investigated  a  sequence  of  configurations  which  would  lead  us  from  one 
equilibrium  configuration  to  another.  The  principal  atoms  driving  the  migration  process 
were  moved  in  small  steps  of  order  0.T0.2  A  with  the  application  of  suitable  constraints, 
while  all  the  other  atoms  were  allowed  to  relax  to  equilibrium. 

We  were  guided  in  our  choice  of  pathways  by  the  recent  comprehensive  investigation  of 
the  diffusion  of  an  0  atom  in  Si.^®  Here,  we  summarize  the  results  of  that  investigation, 
to  motivate  the  bond-breaking  and  bond-making  processes  described  in  the  subsequent 
discussion.  The  migration  of  the  0  atom  is  an  extremely  complex  process  with  both  the 
0  atom  and  a  Si  neighbor  being  required  to  surmount  large  energy  barriers.  The  process 
occurs  through  a  multiplicity  of  paths,  all  of  which  have  an  activation  barrier  of  about  2.5 
eV.  One  path  is  shown  schematically  in  Fig.  1  where  three  Si  atoms  are  labelled  Si(L),  Si(C) 
and  Si(R)  respectively.  The  0  atom  moves  from  the  bond-center  position  between  Si(L)  and 
Si(C)  to  the  bond-center  position  between  Si(R)  and  Si(C).  The  calculated  2.5  eV  activation 
energy  of  migration  for  the  process  is  made  up  of  two  major  contributions: 

(i)  ~  0.6  eV  to  make  the  O  atom  threefold  coordinated  in  an  intermediate  configuration, 
as  shown  in  Fig.  1  (b); 

(ii)  ~  1.6  eV  required  to  allow  Si(C)  to  break  its  bond  . with  Si(R)  and  re-form  its  bond 
with  Si(L),  which  requires  the  system  to  go  through  an  intermediate  configuration  where 
Si(C)  is  only  three-fold  coordinated  (with  a  dangling  bond),  as  shown  in  Fig.  1  (c). 

We  now  describe  a  pathway  for  self-enhanced  G  diffusion,  through  the  migration  of  the 
0  dimer.  This  is  shown  schematically  in  Fig.  2.  The  principal  atoms  driving  the  migration 
process  are  the  two  0  atoms  and  two  Si  atoms,  Si(C)  and  Si(R),  with  Si(L)  and  Si(S)  playing 
a  secondary  role.  In  the  equilibrium  configuration  depicted  in  Fig.  2(a),  the  two  0  atoms 
in  the  dimer  occupy  neighboring  bond-center  positions  with  a  binding  energy  of  0.3  eV,  in 
agreement  with  recent  theoretical  work.^^  All  four  Si-0  bond  lengths  are  ~  1.65  A.  First, 
both  0  atoms  move  in  concert  so  that  they  become  threefold  coordinated  at  the  same  time. 
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as  shown  in  Fig.  2(b).  Here  the  principal  atoms  lie  in  a  four-membered  ring,  with  each  0 
atom  at  the  nominal  midpoint  of  its  migration  path  (in  the  absence  of  the  other  0  atom), 
vertically  above  either  Si(C)  or  Si(R).^^  The  energy  of  this  configuration  is  ~  1.25  eV  above 
equilibrium.  This  is  nearly  twice  the  energy  required  to  make  the  isolated  0  atom  threefold 

coordinated.  The  strong  Si-0  bonds,  drawn  with  solid  lines,  are  ~  1.7  A  in  length,  while 

\  ■ 

the  strained  ones,  drawn  with  dashed  lines,  are  1.9  A  in  length. 

The  next  event  is  the  gradual  re-formation  of  the  bond  between  Si(C)  and  Si(L),  with 
the  other  principal  atoms  moving  to  keep  the  four-membered  ring  intact. As  Si(C)  moves, 
the  0  atoms  rearrange  themselves  so  that  Si(C)  remains  almost  four-fold  coordinated  with 
one  strong  Si-0  bond  of  length  1.7  A  and  another,  considerably  strained  Si-0  bond  of  length 
2.1  A.  The  bond  between  Si(R)  and  Si(S)  becomes  elongated  during  this  motion  of  Si(C). 
The  configuration  with  the  highest  energy  along  this  pathway,  1.5  eV  above  equilibrium, 
is  shown  schematically  in  Fig.  2(c).  Here  Si(C)  is  letting  go  of  the  second  O  atom  while 
the  bond  with  Si(L)  is  close  to  being  formed.  After  the  bond  between  Si(C)  and  Si(L)  is 
re-formed,  the  forces  on  the  two  0  atoms  are  such  that  when  they  move  freely,  the  system 
spontaneously  rearranges  itself  to  attain  the  final  configuration  shown  in  Fig;  2  (d).  The  low 
1.5  eV  activation  energy  is  due  to  the  presence  of  the  second  0  atom  which  enables  Si(C) 
to  be  almost  fourfold-coordinated  during  the  migration  process. 

We  now  show  how  H-enhanced  0  diffusion  takes  place  through  the  migration  of  oxygen- 
hydrogen  complexes.  We  studied  the  equilibrium  configuration  of  the  oxygen- hydrogen 
complex  in  the  positive,  neutral  and  negative  charge  states.  Three  configurations  of  the 
complex  are  illustrated  in  Fig.  3.  In  Fig.  3(a),  the  0  atom  and  the  H  atom  occupy 
neighboring  bond-center  sites.  We  denote  this  configuration  by  BC.  In  Fig.  3(b),  the  H 
atom  occupies  an  anti  bonding  position  relative  to  a  Si  atom  bonded  to  an  0  atom,  which  we 
denote  by  AB.  In  Fig.  3(c),  the  H  atom  is  directly  bonded  to  the  0  atom,  which  we  denote 
by  OH.  We  find  that  the  stability  of  each  configuration  and  the  nature  of  the  equilibrium 
configuration  is  strongly  dependent  on  the  charge  state. 

•  positive  charge  state:  BC  is  the  equilibrium  structure,  and  is  bound  by  0.5  eV 
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compared  with  isolated,  bond-centred  0  and  H"*"  impurities.  OH  is  about  0.1  eV  higher  in 
energy,  and  is,  therefore,  also  strongly  bound  compared  to  the  same  isolated  impurities.  AB 
is  unstable  with  respect  to  spontaneously  distorting  into  BC. 

•  neutral  state:  BC  is  the  equilibrium  structure,  and  is  bound  by  0.3  eV  compared  with 

isolated,  bond-centred  0  and  impurities.  AB  is  0.2  eV  higher  in  energy  and  marginally 

\ 

bouna,  while  OH  is  about  0.5  eV  higher  in  energy  and  is  unbound,  in  comparison  with  the 
isolated  neutral  mpurities. 

•  negative  charge  state:  AB  is  the  equilibrium  structure,  and  is  bound  by  0.4  eV 
compared  with  isolated,  bond-centred  0  and  the  intersititial  H“  impurity.  OH  and  AB  are 
both  much  higher  in  energy  (over  1  eV)  and  are  thus  unbound  with  respect  to  the  isolated 
impurities. 

The  above  results  have  a  logical  consistency  with  the  published  results  for  H  in  Si.^’^  In 
both  cases,  H+  and  H°  prefer  a  site  midway  between  two  Si  neighbors,  where  the  electronic 
density  is  high.  When  the  H+  and  H°  are  placed  next  to  a  Si-O-Si  bridge,  the  Si-H-Si  unit 
buckles  a  little.  Finally,  H“  in  pure  Si  prefers  the  tetrahedral  interstitial  site,  where  the 
electron  density  is  low.  When  placed  next  to  a  Si-O-Si  bridge  the  H~  moves  towards  one  of 
the  Si  neighbors  of  the  0  atom  to  attain  the  AB  configuration.  The  OH  configuration  is  not 
the  equilibrium  configuration  for  any  charge  state,  probably  because  it  makes  the  O  atom 
three-fold  coordinated. 

A  pathway  for  the  migration  of  the  oxygen-hydrogen  complex  in  the  neutral  state  is  shown 
The  equilibrium  BC  configuration  is  shown  in  Fig.  4  (a).  During  the  migration 
process,  the  two  impurity  atoms  move  in  the  following  sequence.  First,  the  0  atom  moves 
towards  the  Si-Si  bond  along  which  the  H  atom  sits,  to  take  up  the  bond-center  position 
between  Si(C)  and  Si(R).^‘*  Owing  to  the  strain  created  by  the  proximity  of  the  0  atom,  the 
H  atom  moves  into  an  antibonding  position  adjacent  to  Si(C),  as  shown  in  Fig.  3(b).  This 
configuration  has  an  energy  only  0.3  eV  higher  than  equilibrium.  The  maximum  increase  in 
energy  during  this  stage  is  only  0.55  eV.  Next,  Si(C)  moves  to  re-form.its  bond  with  Si(L), 
while  the  H  atom  keeps  adjusting  its  position  to  maintain  a  strong  bond  with  Si(C).  When 
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the  bond  between  Si(L)  and  Si(C)  is  almost  re-formed,  in  order  to  complete  the  migration 
process,  the  H  atom  needs  to  move  into  a  position  in-between  Si(R)  and  Si(S).  As  the  H 
atom  climbs  the  energy  barrier  to  accomplish  this,  it  goes  first  into  an  antibonding  position 
with  Si(R).  The  maiximum  energy  during  the  whole  process  corresponds  to  the  configuration 
shown  in  Fig.  4  (c),  which  has  an  energy  1.2  eV  above  equilibrium.  The  low  1.2  e  V  activation 
energy  for  the  migration  of  the  oxygen-hydrogen  complex  is  a  result  of  two  reasons:  (i)  the 
presence  of  the  H  atom  opens  up  the  bond  between  Si(C)  and  Si(R),  making  it  easy  for  the 
O  atom  to  move  in-between;  (ii)  the  H  atom  keeps  rearranging  itself  to  keep  Si(C)  nearly 
fourfold  coordinated,  throughout  the  migration  process. 

We  now  compare  our  results  with  earlier  theory.  For  the  equilbrium  structure  of 
the  neutral  oxygen-hydrogen  complex,  our  results  are  in  agreement  with  Estreicher^'*  and 
in  disagreement  with  Jones,  et  al.^®.  The  energy  surfax^e  proposed  by  Estreicher  is  of 
questionable  quantitative  accuracy,  given  his  high  estimate  of  4.1  eV  for  the  migration 
energy  of  the  O  atom,  compared  to  more  accurate  estimate  of  of  2.6  eV  of  Jones  et  al.  The 
pathways  proposed  by  b.oth  investigators  for  the  migration  of  the  oxygen-hydrogen  complex 
requires  the  complex  to  break  up  and  the  H  atom  to  diffuse  independently  of  the  0  atom. 
The  pathway  we  propose  shows  that  this  is  unnecessary. 

In  the  case  of  self-enhanced  0  diffusion,  the  pathway  that  we  found  has  some  features 
in  common  with  that  proposed  by  Snyder  et  al.^^  They  obtained  a  four-membered  ring-like 
structure,  like  that  shown  in  Fig.  2  (b),  during  the  migration  of  the  0  dimer.  However, 
their  semi-empirical  method  gave  them  a  completely  incorrect  energy  surface.  They  found 
the  ring-like  configuration  to  be  metastable  and  nearly  identical  in  energy  with  that  of  the 
equilibrium  configuration.  Our  calculations,  however,  reveal  that  the  ring-like  configuration 
is  almost  1.3  eV  higher  in  energy  than  equilibrium  and  is  unstable. 

In  conclusion,  we  have  found  mechanisms  for  enhanced  0  diffusion  in  Si.  Self-enhanced 
0  diffusion  proceeds  via  the  migration  of  0  dimers.  The  calculated  activation  energy  is  1.5 
eV,  which  is  close  to  the  experimental  estimate'*”®  of  1.7  eV  for  the  activation  energy  of 
thermal  donor  formation.  Thus,  we  conclude  that  the  formation  of  thermal  donors  occurs 
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through  the  migration  of  0  dimers.  H-enhanced  0  diffusion  occurs  through  the  migration 
of  the  oxygen-hydrogen  complex.  The  calculated  diffusion  activation  energy  is  ~  1.2  eV, 
again  in  good  agreement  with  experiment. In  each  process,  the  reduction  in  the  activation 
energy  by  about  1.0  eV,  in  comparison  with  the  2.5  eV  activation  energy  of  the  0  atom,  is 

the  result  of  the  saturation  of  a  dangling  bond  on  one  of  the  principal  Si  atoms  during  the 

\ 

migration  process. 
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This  was  achieved  by  moving  Si(C)  towards  Si(L),  while  the  other  three  principal  atoms 
were  still  constrained  as  before. 

During  this  motion  of  the  0  atom,  Si(C)  was  moved  into  a  position  on  the  dotted  line 
parallel  to  the  <  001  >  direction,  as  shown  in  Fig.  4,  and  was  allowed  to  relax  only  in  this 
direction.  All  other  atoms  were  allowed  to  relax  to  equilibrium. 
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FIGURES 

FIG.  1.  The  migration  of  the  0  atom  in  Si.  (a)  initial  equilibrium  configuration,  (b)  0  at  the 
mid-point  of  the  path,  (c)  saddle-point  configuration,  (d)  final  equilibrium  structure.  The  small 
dots  are  the  nominal  sites  of  Si  atoms  in  the  bulk.  The  open  circles  are  the  current  positions  of 
the  left  (L),  central  (C)  and  right  (R)  Si  atoms.  The  large  shaded  circle  is  the  0  atom. 

FIG.  2.  The  migration  of  the  0  dimer,  (a)  initial  equilibrium  configuration,  (b)  intermediate 
four-membered  ring  structure,  (c)  saddle-point  configuration,  (d)  final  equilibrium  structure.  The 
dotted  lines  indicate  constraints.  See  Refs.22  and  23. 

FIG.  3.  The  stable  configurations  of  the  oxygen-hydrogen  complex,  (a)  adjacent  bond-center 
configuration  (BC).  (b)  adjacent  anti-bonding  configuration  (AB).  (c)  OH  configuration. 

FIG.  4.  The  migration  of  the  oxygen-hydrogen  complex  in  the  neutral  charge  state,  (a) 
initial  equilibrium  configuration,  (b)  intermediate  four-membered  ring  structure,  (c)  saddle-point 
configuration,  (d)  final  equilibrium  structure.  The  dotted  line  indicates  constraints.  See  Ref.24. 
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Abstract 

We  report  molecular  dynamics  simulations  of  nonequilibrium  heat  flow  in  a  solid 
system  in  the  local-equilibrium  or  hydrodynamic  approximation,  and  demonstrate  that 
local  equilibrium  can  be  achieved  in  a  surprisingly  small  number  of  atomic  layers.  From 
the  dynamical  simulations,  we  calculate  the  thermal  boundary  (Kapitza)  conductance  that 
arises  from  heat  flow  across  Si  grain  boundaries  and  compare  with  the  traditional 
approach  based  on  calculating  the  transmission  and  reflection  of  hannonic  phonons  at  the 
grain  boundary  interface. 

PACS:  66.70.+f,  44.30.+V,  61.72.Mm,  05.60.+W 
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Molecular  dynamics  has  been  a  powerful  method  to  study  the  dynamical 
evolution  of  systems  of  atoms,  namely  a  gas.  a  liquid,  or  a  soUd.  The  positions  qi  and 
momenta  pi  of  all  the  atoms  in  the  system  are  the  dynamical  variables  whose  evolution  in 
time  is  governed  by  Hamilton's  equations  of  motion.  The  classical  Hamiltonian  (total 
energy)  of  the  system  is  detcnnined  from  a  set  of  interatomic  potentials  (conventional 
molecular  dynamics  1 1))  or  ftom  density-functional  theory  (Car-PanineUo  dynamics  (21). 

In  principle,  simulations  are  done  by  assuming  one  of  Gibbs’  statistical  mechanical 
ensembles  that  specify  which  externally  controlled  macroscopic  variables  are  kepi 
constant.  The  resulting  qi(i)  and  Pi(i)  are  the  trajectories  in  phase-space  of  one  replica  in 
the  ensemble.  Typic<tily.  one  assumes  a  uniform  and  constant  temperature  and  either  a 
constant  volume,  or  fi  constant  pressure  or  stress.  The  tr^ectories  are  then  computed  with 
the  consiramts  imposed  by  the  values  of  the  externally  controlled  variables. 

When  extemslly  controlled  macroscopic  variables  such  as  temperature  or  pressure 
arc  non-uniform  over  a  sample,  the  system  is  in  a  nonequilibrium  state  and  irreversible 
transport  occurs.  For  homogeneous  systems,  such  as  liquids  or  uniform  solids,  linear 
transport  coefficients  can  be  computed  from  equilibrium  simulations  (3-5]  using  the 
Kubo  formalism  (6).  This  formalism,  however,  is  not  useful  for  inhomogeneous  systems 
ie.g.  a  polycrystal,  or  an  interface  between  two  different  solids.).  For  such  systems,  one 
must  perfonn  nonet^uilibrium  simulations.  Such  nonequilibrium  simulations,  however, 
present  different  challenges  that  have  not  been  systematically  addressed  so  far.  For 
example,  for  noneqiiilibrium  heat  flow,  we  ate  aware  of  only  one  attempt  to  establish  a 
steady  state  in  a  honiogeneous  Lennard-Jfones  liquid  (7]. 

In  this  Letter,  we  report  the  first  nonequilibrium  dynamical  simulation  of  heat 
flow  in  a  solid  when  a  temperature  difference  is  externally  maintained  between  two 
planes.  First,  we  sl  iow  that  it  is  possible  to  establish  a  uniform  lempcraiuic  gradient 
between  the- two  planes  in  the  case  of  a  single  crystal.  These  results  establish  that  local 
equilibrium  can  be  achieved  in  a  surprisingly  small  number  of  atomic  layers.  Then,  we 
simulate  heat  flow  in  a  bicrystal  and  demonstrate  that  a  temperature  discontinuity  is 
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established  across  the  grain  boundary.  In  terms  of  this  discontinuity  we  are  able  to 
compute  the  Kapitza  conductance  {8,  9]  of  the  grain  boundary.  We  report  such 
calculations  for  two  different  grain  boundaries  in  Si  using  the  Stilling©r-Weber  classical 
three-body  potential  (10]. 

Traditional  calculations  of  the  Kapitza  conductance  have  been  based  on  the  lattice 
dynamical  theory  in  which  one  computes  the  transmission  coefficient  of  harmonic 
phonons  at  the  interface  [11-14]  (see  eqn.  (6)  below).  For  this,  one  needs  to  match 
wavevectors  and  modes  of  phonons  with  the  same  frequency  on  the  incident  and  the 
transmitted  sides  of  the  interface.  The  number  of  such  matching  conditions  is  governed 
by  the  number  of  atoms  in  the  surface  unit  cell  containing  a  periodic  segment  of  the 
interface  parallel  to  the  boundary  plane.  Although  the  theoretical  formalism  of 
wavevector  matching  for  a  general  interface  exists  [12],  it  has  so  far  been  applied  to 
examples  with  only  a  few  atoms  in  the  surface  unit  cell.  However,  even  for  the  simplest 
case  of  a  twin  boundary,  there  are  many  atoms  in  the  surface  unit  cell,  and  using  the 
traditional  mode  counting  approach  becomes  a  nontrivial  exercise.  The  molecular 
dynamics  method  to  calculate  Kapitza  conductance  is  superior  to  the  traditional  method 
because:  (i)  it  is  much  simpler  to  set  up  for  any  crystal  structure  and  any  interface;  (ii)  the 
accuracy  of  the  calculations  can  be  systematically  increased  by  merely  increasing  the 
supercell  size  and  the  total  simulation  time;  and  (iii)  all  anharmonic  effects  are  included 
automatically. 

Before  describing  the  simulation  details,  it  is  relevant  to  address  the  major 
obstacles  to  a  successful  nonequilibrium  simulation,  that  would  help  explain  why  not 
many  such  simulations  have  been  reported  to  date.  The  major  difficulty  stems  from  the 
fact  that  in  order  to  establish  an  accurate  thermal  profile,  one  needs  to  partition  the 
supercell  into  various  sections  and  define  a  local-equilibrium  temperature  for  each 
section.  To  accurately  define  a  local  temperature,  one  would  like  to  have  sections  that  are 
long  compared  to  the  mean  free  path  (mfp)  [5]  of  the  heat  carriers,  i.e.,  phonons  in  our 
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case.  Starling  from  n  Boltzmann  transport  theory  formalism  (15),  one  obtains  the 
following  expression  for  the  averaged  phonon  mfp  for  an  isotropic  solid; 

p  Cv 

where  00  is  the  bulk  ihennal  cemhictivity.  p  the  mass  density  of  the  crystalline  solid.  C, 
she  phonon  contribution  to  the  speciftc  heat  at  constant  volume,  and  <Vph>  an  average 
phonon  velocity-  The  last  two  quantities  are  defined  by: 

C,  .  Z  J  d3(J  hoKX.  ;  P> 

X 


Z  \ d’q  0(X.  qX 


Zjd’qo(k,q)^^ 


X 

where  (i)(X,  q)  is  the  phonon  frequency  corresponding  to  the  mode  X  and  wavevcctor  q, 
f(<i),  T)  is  the  pbonoii  occupation  given  by  the  Bose-Einstein  factor,  and  all  integrations 
are  defined  over  the  whole  Brillouin  Zone.  We  have  calculated  the  phonon  spectrum  of  Si 
using  the  Stillinger-Weber  potential  [10)  and  obtained  values  of  Cv  and  <Vph>.  The 
resulting  <lph>  is  8h.>wn  in  Pig.  1  as  a  function  of  temperature.  It  is  clear  that  the  phonon 
mfp  is  large,  from  20  ran  to  several  microns.  Even  at  the  Debye  temperature  (0d~  ^50  K 
(16))  the  -  30  ran.  Thus,  supercells  that  arc  large  compared  with  the  phonon  mfy 
should  contain  thousands  of  atoms.  Another  measure  of  the  difficulty  is  the  average 
phonon  lifetime  Xp^  -  <lph>/<Vph>v  which  is  -  0.01  ns  at  T  -  and  orders  of 
magnitude  larger  at  lower  temperatures.  Again,  in  order  to  generate  a  good  steady-state 
profile,  one  would  expect  simulations  to  be  longer  than  tph  which  might  be  prohibitively 

long  at  low  temperatures. 
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In  spite  of  great  strides  in  the  understanding  of  nonequilibrium  statistical 
mechanics  [17],  there  is  no  established  theory  to  systematically  predict  the  smallest 
supercell  size  and  the  shortest  simulation  time  required  to  establish  a  thermal  profile  with 
reasonable  accuracy.  From  a  practical  point  of  view,  it  is  therefore  important  to  explore 
these  questions  with  computer  experiments  on  supercells  of  various  sizes  and  various 
simulation  times.  From  the  foregoing  discussion  it  would  seem  that  simulations  in  Si  are 
completely  unfeasible  at  very  low  temperatures  (T  «  ©d)-  I*  seems  to  be  a  significantly 
daunting  task  even  at  T  ~  ©d-  As  a  first  step  towards  a  more  general  investigation  we 
have,  therefore,  concentrated  only  at  a  temperature  range  just  below  ©d*  Results  of  three 
such  simulations  reported  below  are  quite  encouraging:  it  is  possible  to  establish  steady- 
state  thermal  profiles  using  sections  that  are  much  smaller  than  the  phonon  mfp  with 
reasonably  long  simulation  times. 

We  now  describe  details  of  the  simulation  set  up.  All  supercells  are  in  the  shape 
of  rectangular  parallelopipeds  (Fig.  2)  with  a  length  of  a  few  tens  of  nm  in  the  direction  of 
heat  flow,  taken  normal  to  the  grain  boundary  plane  (x-axis  of  Fig.  2).  A  few  A  of 
vacuum  is  left  beyond  the  two  end  faces  parallel  to  the  gram  boundary  plane;  otherwise 
the  application  of  a  temperature  gradient  would  be  incompatible  with  periodic  boundary 
conditions  employed  in  the  simulation.  The  y  and  z  dimensions  of  the  supercell  are  set 
equal  to  the  dimensions  of  one  periodic  segment  of  the  grain  boundary  plane.  For 
simulations  on  (i)  uniform  Si,  (ii)  the  Z-5  twin  boundary,  and  (iii)  the  2=:13  twin 
boundary,  these  dimensions  are  a  x  a,  V5/2  a  x  a  and  Vl3/2  a  x  a  respectively,  a  being  the 
Si  lattice  constant.  For  simulations  with  the  twin  boundaries,  the  two  grains  are  chosen  of 
equal  length  so  that  the  grain  boundary  is  located  precisely  in  the  middle  of  the 
simulation  cell.  In  order  to  define  a  temperature  profile,  the  sample  solid  is  equally 
divided  into  an  odd  number  of  sections  parallel  to  the  grain  boundary  plane,  with  the 
grain  boundary  fully  contained  in  the  middle  section.  For  each  section  S,  an  average 
temperature  T(S)  is  defined  by: 
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(4) 


where  Ns  is  the  numher  of  particles  contained  in  the  section  S,  N  the  total  number  of 
particles  in  the  system,  and  o  denotes  averaging  over  the  simulation  time.  The  quantity 
on  the  l.h.s.  of  cq.  (4)  is  the  time-averaged  total  kinetic  energy  of  section  S,  while  the  sum 
on  the  r.h.s.  of  eq.  (4)  is  the  total  energy  of  the  phonons  corresponding  to  an  equilibrium 
temperature  T(S).  TTie  zero-point  energy  has  been  excluded  from  the  phonon  energy 
expression,  so  that  T(S) «  0  corresponds  to  static  atoms.  In  the  high  temperature  limit,  i.e. 
for  T(S)  »  Od,  r  h-s-  of  3Ns1cbT(S)/2,  the  usual  definition  of  temperature 

employed  in  most  simulations.  For  consistency,  wc  have  used  eq.  (4)  with  the  phonon 
dispersion  aCX.  q)  olitained  from  the  StiUinger-Weber  potential  (10). 

The  thermal  gradient  is  applied  by  maintaining  the  two  end  sections  at  constant 
but  different  tempwMures  Tj  and  T2-  For  isothermal  simulations  of  the  two  end  sections, 
we  employ  Hoover’s  method  [18J  in  which  a  velocity-dependent  drag  term  is  introduced 
into  the  equation  of  motion.  This  method  of  constant  temperature  is  more  suitable  for  the 
present  simulations  than  the  canonical  method  of  No6  [19]  in  which  there  is  an  inherent 
temperature  fluctuation  of  order  T/VNi,  Nj  being  the  number  of  particles  in  section  S,  and 
T  the  desired  temperature.  In  all  sections  other  than  the  two  ends,  the  atoms  move 
according  to  Ncwtoii's  equation  of  motion  with  forces  derived  from  the  Stillinger-Weber 
potential.  The  equa  tions  of  motion  are  numerically  integrated  by  the  Leapfrog  Verlci 
algorithm  [1]  using  a  time-step  of  0.8  femtoseconds. 

Fig.  3  displays  ±t  tiiennal  profile  for  crystalline  Si  (i.«.  no  grain  boundary)  after  a 
total  simulation  time  of  1.0  nanosecond,  which  was  found  to  be  adequate  to  establish  a 
steady-state  profile  (sec  discussion  below).  The  supercell  is  of  length  45  nm,  and  the 
sample  is  divided  into  21  equal  sections  each  containing  30  atoms.  The  yz  plane  is  taken 
to  be  the  [100]  plane  of  the  Si  crystal.  The  cold  and  the  hot  ends  Le.  section  #s  1  and  21 
are  maintained  at  constant  temperatures  of  520  K  and  630  K  respectively  throughout  the 
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entire  simulation.  Even  this  unrealistically  large  temperature  gradient  yields  a  linearly 
varying  profile  across  a  large  part  of  the  sample  (section  #s  3  through  16),  in  accordance 
with  Fourier’s  law.  The  end  effects,  also  present  in  a  previous  simulation  for  Lennard- 
Jones  liquid  [10],  are  due  to  enhanced  scattering  of  phonons  arising  from  the  additional 
"drag"  terms  introduced  to  maintain  constant  temperatures  [18]  at  the  end  sections.  The 
asymmetry  about  the  middle  section  (#  11)  is  due  to  a  decrease  in  the  thermal 
conductivity  of  Si  with  increase  in  temperature[20, 21]  in  the  range  considered. 

Rg.  4  displays  a  similar  thermal  profile  when  the  middle  section  (#  1 1)  contains  a 
symmetric  2-13  tilt  boundary,  with  a  supercell  of  length  25  nm,  containing  840  atoms 
(the  profile  for  a  2-5  boundary  is  qualitatively  very  similar,  and  is  not  shown  explicitly). 

This  profile  is  dramatically  different  from  the  pure  Si  case  (Fig.  3)  in  that  a  sharp 
discontinuity  in  temperature,  (AT)gb  appears  across  the  grain  boundary.  In  addition,  there 
axe  pronounced  end  effects  due  to  the  same  reasons  as  for  the  pure  Si  case.  However,  the 
sample  was  chosen  large  enough  to  extract  the  linear  slope  (dT/dx)xtal  in  the  crystalline 
material  away  from  the  ends,  as  given  by  the  slopes  of  the  two  fitted  lines  AB  and  CD. 

The  Kapitza  conductance  can  be  readily  calculated  in  terms  of  the  bulk  thermal 
conductivity  gq  using  the  formula: 

(5) 

(aT)gb 

The  values  obtained  by  eq.  (5)  from  the  lower  (AB)  and  upper  (CD)  slopes  (in 
conjunction  with  oq  values  for  the  corresponding  lower  and  higher  temperatures)  are 
averaged  to  yield  gr  “  0-8  GW/(m?K)  for  the  2=13  boundary.  A  similar  calculation  for 
the  2=5  boundary  yields  Gk  =  0.9  GW/(m^K). 

The  smooth  thermal  profiles  of  Figs.  3  and  4  are  surprising  considering  the  fact 
that  each  section  is  only  1-2  nm  long,  whereas  the  mJ^  ~  30  nm!  In  order  to  obtain  a 
deeper  insight  into  why  it  works,  it  is  necessary  to  examine  the  profile  as  it  develops  with 
time.  Monitoring  of  the  time  evolution  of  the  profiles  reveals  that  signatures  of  high 
temperature  gradients,  at  the  ends  and  across  the  grain  boundary,  get  established  much 
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earlier  than  the  flatter  profiles  in  the  bulk  '’ciystalline''  regions.  For  the  simulation  on  the 
51=13  boundary.  Fig.  5  shows  die  tempemure  of  a  given  section  (#  7,  chosen  for 
concreteness)  as  a  furicdon  of  the  simulation  time.  Rg.  5  also  displays  the  time  evolution 
of  the  temperature  discontinuity  (AT)gB  across  the  grain  boundary.  Both  quanUties 
oscillate  significantly  for  the  first  few  hundred  ps.  and  attain  steady-state  values  in  amc 
of  order  a  ns.  Howe\er,  HMt  fractional  fluctuation  of  (AT)gb  about  its  steady-state  value 
reduces  to  less  than  10%  only  after  70  ps.  It  takes  much  longer,  at  least  a  ns,  to  achieve 
the  same  order  of  nccuracy  in  the  slope  (dT/dx)xui,  because  it  entails  an  absolute 
fluctuation  in  T(S)  less  than  a  fraction  of  1  K,  for  all  S  in  the  crystalline  region.  The 
observation  that  thermal  profiles  get  established  faster  in  regions  of  larger  phonon 
scattering  rates,  implies  that  local  thermal  equilibrium  is  established  presumably  by  a 
certain  critical  numlcr  of  phonon  scaiiering  events,  even  if  the  section  is  much  smaller 
than  the  mjp.  Using  the  Boltemann  relaxation  approximation,  and  an  energy  independent 
relaxation  time,  we  have  estimated  the  number  of  phonon  scattering  events  witiiin  each 
section  (in  the  ciystilline  region)  to  be  a  few  thousand  for  a  total  simulation  time  of  1  ns. 
This  number  seems  to  be  suffident  to  establish  a  smooth  thermal  profile. 

In  order  to  recast  the  above  values  of  Kapitza  conductance  in  terms  of  a  lattice 
dynamical  model,  we  recall  that  or  can  be  expressed  as  the  temperature  derivative  of  the 
phonon  heat  ciurent  density  across  flic  interface  (6): 

I  q)  Vx(X,  q) 

vx>0 

where  i(?i,  q)  is  the  transmission  coefficient  of  the  (X.  q)  phonon,  i.e.  the  fraction  of  the 
energy  of  an  incident  phonon  of  wavevector  q  and  mode  X  that  gets  transmitted  through 
the  boundary,  and  vx  is  the  phonon  group  velocity  normal  to  the  boundary.  The 
integmdon,  as  befoie,  is  over  the  entire  Bfillouin  Zone,  and  the  condition  Vx  >  0  chooses 
only  those  phonow  that  are  incident  on  the  grain  boundary  from  the  left  side.  As 
indicated  earlier,  calculating  t(X,  q)  from  wavevector  matching  equations  is  not 
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Straightforward  for  the  case  of  a  grain  boundary.  However,  once  or  values  are  obtained 
by  MD  simulations,  one  could  estimate  an  average  transmission  coefficient  <t>  =  ok  / 1, 
where: 

/ 


i-S 

'  X 


vx>0 


q)  Vx(X,  q) 


af(G),T) 

3T 


(7) 


The  quantity  I  can  be  interpreted  as  the  maximum  possible  value  of  ok.  which  occurs 
under  perfect  phonon  transmission  across  the  interface,  i.e.  t(X,  q)  =  1  for  all  X  and  q. 
Using  the  Stillinger-Weber  phonon  dispersion  at  our  mean  simulation  temperature  of  575 
K  we  obtain  I  =  1.4  GW/Cm^K),  which  yields  <t>  =  0.65  and  0.57  for  the  S-5  and  the 
Z=13  boundaries,  respectively.  The  fact  that  <t>  <  1  for  both  the  boundaries  is  a 
consistency  check  on  the  MD  simulations.  The  relatively  large  values  of  the  average 
transmission  coefficient  are  consistent  with  the  fact  that  all  Si  atoms  in  the  two  grain 
boundaries  are  fourfold  coordinated  [22, 23],  and  therefore  the  local  phonon  modes  in  the 
Interface  region  are  not  expected  to  be  very  different  from  that  in  crystalline  Si.  A  more 
systematic  study  of  the  variation  of  Ok  with  the  grain  boundary  tilt  angle  is  underway. 

In  summary,  we  have  demonstrated  that  local  equilibrium  can  be  achieved  in  a 
nonequilibrium  molecular  dynamics  simulation  on  sections  that  are  much  smaller  than  the 
phonon  m^.  The  technique  is  illustrated  by  establishing  steady-state  thermal  profiles  in 
crystalline  Si  and  across  two  different  symmetric  tilt  grain  boundaries  in  Si.  The  results 
establish  the  validity  of  Fourier’s  law  even  under  very  large  temperature  gradients  (of 
order  1  K/A),  and  give  rise  to  dramatic  temperature  discontinuities  across  interfaces 
between  two  different  crystals.  This  allows  for  the  first  time  a  simple  method  to 
determine  the  thermal  boundary  (Kapitza)  conductance  associated  with  a  general  grain 
boundary,  and  opens  the  door  for  more  nonequilibrium  simulations  in  solid  and  liquid 
systems. 

This  research  was  supported  in  part  by  Lockheed  Martin  Energy  Research  Corp. 
under  DOE  contraa  No.  DE-AC05-96OR22464,  and  ONR  grant  No.  NOOO 14-95- 1-0906. 
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Figure  capiioM 


1 .  Avenge  plwnon  mean  fiee  pa*  <>i*>  «  crystalline  Si  as  a  fiinction  of  tempetamie. 
See  text. 

2.  Schematic  she  lving  the  shape  of  supcrcell  employed  in  all  simulations.  The 
superccU  is  divided  into  an  odd  number  (13  in  this  figure)  of  equal  sections  parallel 
to  the  yz  plane.  The  grain  boundary  is  entirely  contained  in  the  middle  section 
(section  7).  Periodic  boundary  condition  is  employed  on  all  sides.  Some  vacuum  is 
left  beyond  the  two  end  sections  which  are  kept  at  two  different  temperatures  Ti 

andT2. 

3.  Thermal  piofib  for  pure  Si.  The  superoell  is  450  A  along  the  x-axis  and  contains 
630  atoms  (21  sections,  each  containing  30  atoms).  The  IlOOJ  plane  of  the  Si  crystal 
is  oriented  parsllel  to  the  yz  plane.  Total  simulation  time  is  1  ns. 

4.  Thermal  profile  for  a  1=13  twin  boundary.  The  supercell  is  250  A  long  and 
contains  840  atoms.  The  grain  boundary  is  obtained  by  rotating  the  two  "grains"  of 
the  original  crystal  in  opposite  senses  through  an.angle  of  11.31*  about  the  z-axis. 
Total  simulation  time  is  1  ns. 

5.  Time  evolulio  a  during  simulation  on  the  X°13  twin  boundary:  lower  (soUd)  curve 
monitors  the  iMnperatuie  of  section  #7  measured  relative  to  the  cold  end  (520  K); 
upper  (dashed)  line  describes  the  behavior  of  the  temperature  discontinuity  across 

grain  boundary  (rescaled  by  a  factor  of  half). 
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Ihc  tilled  planar  surface  have  also  been  predicted  [10,1 1].  Here,  we  focus 
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